(c)Ferdi Hellweger, Northeastern University, Boston.
NOTE: The applet requires Java 5 or higher. Java must be enabled in your browser settings. Mac users must have Mac OS X 10.4 or higher. Windows and Linux users may obtain the latest Java from Sun's Java site.
QUICK START: Click setup, then go.
powered by NetLogo
view/download model file: DiffusionBias.nlogo
[extract from Hellweger and Bucci, 2009]
Microbes can be motile. Chemotaxis by bacteria using flagella and diel vertical migration by phytoplankton using buoyancy regulation are two examples. These active transport behaviors are straightforward to implement in an IBM. However, microbes are also often transported passively by the ambient hydrodynamics. Advection is simulated by simply applying the velocity to the individual’s position, and diffusion can be simulated using a random walk scheme. For a constant diffusivity, the step size is drawn from a normal distribution with standard deviation based on the diffusivity:
where z (m) is the position (e.g. depth), R is a random number from a standard normal distribution, K (m2 s-1) is the diffusivity and t (s) is time. However, the diffusivity often varies in space. Mixing in the ocean is driven by wind stress at the water surface and the diffusivity generally decreases with depth. In the metalimnion (thermocline) of a lake, density stratification has a stabilizing effect that leads to lower diffusivity. Application of the simple random walk model to this situation leads to an unrealistic net transport towards areas of lower diffusivity. This can be explained theoretically and numerically, or conceptually by the fact that “over a given time interval, particles at a particular location are influenced by slightly more energetic eddies originating in the area of high diffusion compared with somewhat less energetic eddies from the low diffusion region” (Visser, 1997). The problem can be avoided by transforming the coordinate system so that the diffusivity is constant or by using a different formulation for the step size (Hunter et al., 1993; Visser, 1997; Nagai et al., 2003):
where K'(z) is the derivative of K(z). This formulation is different from the simple random walk in that it includes a corrective advection term and the diffusivity is taken at an offset location. Nagai et al. (2003) found that the correction scheme does not work for zero diffusivity and solved this problem by incorporating background diffusivity. Nopass boundaries (e.g. water surface) are typically simulated as reflecting.
1. Adjust the slider parameters (see below), or use the default settings.
2. Press the "setup" button.
3. Press the "go" button to begin the simulation.
4. View the bacteria population develop over time.
5. Turn on the "Correction" and observet the effect.
Try adjusting the parameters. Use a different pattern for the diffusion. Can the correction remove the bias for all patterns and parameters?
Diffusion Bias References:
Hellweger, F. L., Bucci, V. 2009. A bunch of tiny individuals – Individual-based modeling for microbes (review paper). Ecological Modelling, 220(1):8-22.
Hunter, J.R., Craig, P.D., Phillips, H.E., 1993. On the use of random walk models with spatially variable diffusivity. J. Comput. Phys. 106, 366–376.
Visser, A.W., 1997. Using random walk models to simulate the vertical distribution of particles in a turbulent water column. Mar. Ecol. Prog. Ser. 158, 275–281.
Nagai, T., Yamazaki, H., Kamykowski, D., 2003. A Lagrangian photoresponse model coupled with 2nd-orderclosure. Mar. Ecol. Prog. Ser. 265, 17–30.
Wilensky, U. & Reisman, K. (1999). Connected Science: Learning Biology through Constructing and Testing Computational Theories -- an Embodied Modeling Approach. International Journal of Complex Systems, M. 234, pp. 1 - 12. (This model is a slightly extended version of the model described in the paper.)
Wilensky, U. & Reisman, K. (in press). Thinking like a Wolf, a Sheep or a Firefly: Learning Biology through Constructing and Testing Computational Theories -- an Embodied Modeling Approach. Cognition & Instruction.
; diffusionbias.nlogo ; ferdi hellweger ; firstname.lastname@example.org ; 1/12/2010 globals [ dt time dtp ntp xy xy_u xy_i xy_p xy_p_u xy_p_i y_p y_p_u y_p_i y_p_t xt yt x2 y2 dx2 dy2 y3 y4 xp z z_u a_p_u v_p_u maxx minx maxy miny pcolmax pcolmin pcolmid pcolslope DmaxPrev DminPrev PatternPrev D_p_y D_p_t D_p_u D_p_a D_p_b Dp_p_u D_b D_b_u Dp_b Dp_b_u D_b3 D_b_u3 nb_u bn bd bn_a bd_a bd_a_p AnalyticTest mp ] breed [ bguys bguy ] patches-own [ D_p Dp_p bn_a_p ] bguys-own [ ] to setup clear-all set time 0 ;d set dt 0.01 ;d set dtp 0.01 ;d set ntp dt ;d set xy 10 ;cm set xy_u ( xy / 100 ) ;m set xy_i ( max-pxcor - min-pxcor + 1) ;i set xy_p ( xy / xy_i ) ; cm set xy_p_u ( xy_p / 100 ) ;m set xy_p_i 1 ;i set z 1 ;cm set z_u ( z / 100 ) ;m set a_p_u ( xy_p_u * z_u ) ;m2 set v_p_u ( a_p_u * xy_p_u) ;m3 set xp 1e-9 set maxx ( max-pxcor + 0.5 - xp ) set minx ( min-pxcor - 0.5 + xp ) set maxy ( max-pycor + 0.5 - xp ) set miny ( min-pycor - 0.5 + xp ) set pcolmax 17 set pcolmin 12 set pcolmid ( ( pcolmax + pcolmin ) / 2 ) set pcolslope ( ( pcolmax - pcolmin ) / ( Dmax - Dmin ) ) ;random-seed 1 set nb_u nb set-default-shape bguys "circle" if Initial = "Row" [ set yt random-pycor ] if Initial = "Point" [ set yt random-ycor set xt random-xcor ] create-bguys nb_u [ set color black set size 0.5 if Initial = "Spread" [ setxy random-xcor random-ycor ] if Initial = "Row" [ setxy random-xcor yt ] if Initial = "Point" [ setxy xt yt ] ] setD set AnalyticTest false end to setD if not ( Dmax = DmaxPrev ) or not ( Dmin = DminPrev ) or not ( Pattern = PatternPrev ) [ set-current-plot "Diffusion" set-current-plot-pen "D" plot-pen-reset set D_p_y  set y_p  set y_p_i n-values ( max-pycor + 1 ) [?] foreach y_p_i [ set y_p_t ( ? * xy_p ) if Pattern = "Uniform" [ set D_p_t ( Dmax - Dmin ) / 2 ] if Pattern = "Linear" [ set D_p_t Dmin + y_p_t * ( Dmax - Dmin ) / ( xy - xy_p ) ] if Pattern = "Exponential" [ set D_p_t ( 10 ^ ( ( log Dmin 10 ) + y_p_t * ( ( log Dmax 10 ) - ( log Dmin 10 ) ) / ( xy - xy_p ) ) ) ] if Pattern = "Step" [ ifelse y_p_t < ( xy / 2 ) [ set D_p_t Dmin ] [ set D_p_t Dmax ] ] ask patches with [ pycor = ? ] [ set D_p D_p_t ] set D_p_y lput D_p_t D_p_y set y_p lput y_p_t y_p plotxy D_p_t y_p_t ] ask patches [ ifelse pycor = max-pycor [ set D_p_a D_p ] [ set D_p_a [ D_p ] of patch pxcor ( pycor + 1 ) ] ifelse pycor = min-pycor [ set D_p_b D_p ] [ set D_p_b [ D_p ] of patch pxcor ( pycor - 1 ) ] set Dp_p ( D_p_a - D_p_b ) / ( 2 * xy_p ) ] ask patches [ set pcolor ( pcolmax - D_p * pcolslope ) ] set DmaxPrev Dmax set DminPrev Dmin set PatternPrev Pattern ] end to go setD ask bguys [ set D_b [ D_p ] of patch-here set D_b_u ( D_b * 1.0e-5 / ( 100 ^ 2 ) * 86400 ) set Dp_b [ Dp_p ] of patch-here set Dp_b_u ( Dp_b * 1.0e-5 / ( 100 ^ 2 ) * 86400 * 100) set dx2 ( ( random-normal 0.0 1.0 ) * sqrt( 2 * D_b_u * dt ) ) ifelse Correction [ set y3 ( ycor + 0.5 * Dp_b_u * dt / xy_p_u ) if y3 > maxy [ set y3 maxy ] if y3 < miny [ set y3 miny ] set y4 ycor set ycor y3 set D_b3 [ D_p ] of patch-here set ycor y4 set D_b_u3 ( D_b3 * 1.0e-5 / ( 100 ^ 2 ) * 86400 ) set dy2 ( Dp_b_u * dt + ( random-normal 0.0 1.0 ) * sqrt( 2 * D_b_u3 * dt ) ) ] [ set dy2 ( ( random-normal 0.0 1.0 ) * sqrt( 2 * D_b_u * dt ) ) ] set x2 ( xcor + dx2 / xy_p_u ) if x2 > maxx [ set x2 ( maxx - ( x2 - maxx ) ) ] if x2 < minx [ set x2 ( minx - x2 ) ] set y2 ( ycor + dy2 / xy_p_u ) if y2 > maxy [ set y2 ( maxy - ( y2 - maxy ) ) ] if y2 < miny [ set y2 ( miny - y2 ) ] set xcor x2 set ycor y2 ] if time >= ntp [ set-current-plot "Bacteria" set-current-plot-pen "B" plot-pen-reset foreach y_p_i [ set y_p_t ( ? * xy_p ) set bn count bguys-on patches with [ pycor = ? ] set bd ( bn / xy_p ) ;no/cm plotxy bd y_p_t ] if AnalyticTest and Initial = "Row" and Pattern = "Uniform" [ set mp ( nb_u / xy_i / a_p_u ) ;no/m2 ask patches [ set D_p_u ( D_p * 1.0e-5 / ( 100 ^ 2 ) * 86400 ) set y_p_u ( ( pycor - yt ) * xy_p_u ) set bd_a_p ( mp / ( 2 * sqrt ( pi * D_p_u * time ) * exp ( - y_p_u ^ 2 / ( 4 * D_p_u * time ) ) ) ) ;no/m3 set bn_a_p ( bd_a_p * v_p_u ) ] set-current-plot-pen "B_a" plot-pen-reset foreach y_p_i [ set y_p_t ( ? * xy_p ) set bn_a ( sum [ bn_a_p ] of patches with [ pycor = ? ] ) set bd_a ( bn_a / xy_p ) ;no/cm plotxy bd_a y_p_t ] ] set ntp ( time + dtp ) ] set time ( time + dt ) end