Can anyone solve a stochastic differential equation - related to neuroscience research? I'm a neuroscience grad student, and I'm hoping one of ya'll could help me solve this problem regarding particle diffusion. It relates to my research on molecular-level neural plasticity, but I've framed the problem in general terms below. This will probably require you to have some experience with stochastic differential equations. Any insight is much appreciated. Thanks!
The problem in a nutshell

I am ultimately interested in the average steady-state particle density that can be expected in a defined region of space, when the particles are diffusing with random motions -- assume particle diffusion is Brownian motion, with zero drift, and have a constant diffusion rate coefficient.
For example I'm looking for a numerical solution for the average steady-state particle density in two circular subregions with a diameter of 0.8 arbitrary units (au) contained inside a rectangular 3x6 au flat surface (these subregions are equal distance from the edges - see below). Here, let's define 1 au of distance as 1 micron (µm). In this 18 µm² space the particle diffusion coefficient is 0.1 µm²/s, except in the circular subregions; in subregion-1 the diffusion rate is 0.05 µm²/s, in subregion-2 the diffusion rate is 0.01 µm²/s. The boundary conditions can be considered 'rebound', and there are 200 total particles in this closed system. The solution should be able to generalize to any number of particles, but solving this particular situation would be great for starters. 
See the question visually represented here in Fig 1 along with two graphs described below.
Each data-point in the graphs represent the average steady-state particle density for 10 independent simulation trials. 
Monte Carlo Simulation Results

Figures

(Fig 1) I've simulated the above scenario using Matlab. Particles diffused along the 2D surface of a closed rectangular environment that contained two equal-area subregions with a diffusion rate slightly lower than the global diffusion rate of 0.1 µm²/s. The diffusion coefficient (Dcoeff) of the bottom circular subregion was set to 0.01 µm²/s while the top subregion was set to 0.05 µm²/s. There were 500 total time-steps in each individual simulation. The heat-map seen on the right highlights regions of relatively high particle density at the end of the simulation (aka at steady-state); which was done using matrix convolution of particle locations with a Gaussian-shaped mask (top right).


*

*Matlab code for this simulation


(Fig 2) The diffusion coefficients mentioned above resulted in average steady-state particle density values that were vastly different between the two subregions, with little variance (CI envelopes reflect noise for 10 iterations). The top subregion averaged ~15 particles at steady-state (green), while to bottom region averaged ~60 particles (red). The third line (blue) represents the steady-state particle density in a circular region the same size as the two subregions, but was set to the global diffusion rate.


(Fig 3)
This figure shows the effects of holding the Dcoeff of one subregion constant at 0.01 µm²/s while changing the other Dcoeff from 0.01 µm²/s to 0.05 µm²/s (at 0.01 µm²/s increments). 
(left panel) The Dcoeff ratio at each step (calculated by simply dividing the pre-set diffusion rates of the two subregions) and the resulting steady-state particle density in each subregion are almost exactly proportional. 
(right panel) Interestingly, particle availability minimally affects this outcome - as the upper subregion transitioned from 0.01 µm²/s to 0.05 µm²/s it ultimately made an additional 25-30 particles available globally, but only a small portion of them accumulated in the lower subregion (which had a stable Dcoeff of 0.01 µm²/s throughout this ancillary simulation).
The math question and proof of concept

I believe that the particle diffusion can be summarized as Wiener process / Brownian motion with zero drift, and the Kolmogorov forward equation (aka Fokker Plank equation) will describe the time-evolution of the PDF for a random process. However, it's not immediately clear to me how to define the sODE (or sPDE); consequently, I don't even have a good sense of the difficulty-level of this question. Ultimately what I'm looking for, is an equation where I can enter these constants:


*

*Co (Outer box xy dimensions)

*Ci (Inner circle radius)

*Do (diffusion coefficient for Co)

*Di (diffusion coefficient for Ci)

*N  (total number of particles in the closed system)


And the output will provide the expected density of particles in Co and Ci at steady-state. Note that I'm also interested in the rate of change in these values if suddenly D changes from 0.01 µm²/s to 0.05 µm²/s, but first things first.
I believe the steady-state diffusion equation will be of use:
∇D(r)∇n(r)=0
with spatially dependent diffusion coefficient 
D(r)
and prescribed total particle number 
N=∫n(r)dr
If you have an answer, please use an example showing the input/output. I would like to test it against the Monte Carlo simulation to be convinced. If the fact that the two-dimensional geometry has no special symmetry (which prevents a closed-form solution), feel free to just cut the system in half, such that we are only considering the bottom 3x3 square with the circular subregion directly in the center. In fact, if it makes it significantly easier, the outer container can be circular as well. 
Again, any help is much appreciated! Cheers for sharing your skills. 
 A: Alvarez-Ramirez and Meraz have written a paper entitled Assymetric diffusion in heterogenous media. I'm not sure if everybody can access the link, I had to log into my university network first.
In the paper, the following example is considered. Let the domain be $[-L,L]$, i.e. 1D. There are two domains, from $-L$ to $0$ and from $0$ to $L$ with an interface at $x=0$. The paper starts with a descripton of the Monte Carlo simulation. Then the diffusion equations are explained. Consider the diffusivity parameters $D_-$ and $D_+$ for the left and right domain, respectively. Then we have the following equations for the diffusion process:
\begin{align}
\partial_t c(x,t) =D_- \partial_{xx} c(x,t)\quad x\in (-L,0)
\\
\partial_t c(x,t) =D_+ \partial_{xx} c(x,t)\quad x\in (0,L)
\end{align}
with $c(t,x)$ being the concentration. For the interface, I quote: 

Ovaskainen and Cornell provided a rigorous proof that normal diffusion
  across interfaces should meet the following jump condition:

\begin{align}
D_-c(0^-,t) =D_+c(0^+,t)
\end{align}
This formalism can be extend to your scenario easily. I'm currently working on a small implementation for your scenario and will update my answer, once I have produced some results.
