DiffusionBias: Bias in Lagrangian Models with Spatially Variable Diffusion (v1.0)

GETTING STARTED

(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.

MODEL

powered by NetLogo

view/download model file: DiffusionBias.nlogo

WHAT IS IT?

[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.


HOW TO USE IT

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.


THINGS TO TRY

Try adjusting the parameters. Use a different pattern for the diffusion. Can the correction remove the bias for all patterns and parameters?


CREDITS AND REFERENCES

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.

NetLogo References:

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.


PROCEDURES

; diffusionbias.nlogo
; ferdi hellweger
; ferdi@coe.neu.edu
; 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