PDE -- Account for Boundary Conditions

[dhermes]

See: http://www3.ul.ie/gleesonj/Papers/SDEs/AM4062_notes6.pdf

We prefer a " Reflecting boundary condition".

(I added this because I don't want to lose track of the link.)

[bradmonk]

The problem in a nutshell:

We are ultimately interested in the _average steady-state number of particles_ that can be expected in a defined region of space, when the particles are diffusing with random motions -- we assume particle diffusion is Brownian motion, with zero drift, and has a constant diffusion rate coefficient.

For example we want a numerical solution for the average steady-state number of particles 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. Here's the question visually represented...

Monte Carlo Simulation Results

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 can be found here:

• script with syntax highlighting
• script on Github Gist.

The figure below is showing the averages from 10 independent simulation trials. These results are described in a bit more detail below.

(graph above) The diffusion coefficients mentioned above resulted in average steady-state particle number 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 (red), while to bottom region averaged ~60 particles (green). The third line (blue) represents the steady-state particle number in circular region the same size as the two subregions, but was set to the global diffusion rate.

(graph below) 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 number of particles in each synapse 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).