ABMvsLSM: Agent-Based Model vs. Lumped-System Model (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

WHAT'S GOING ON?

Present biogeochemical models typically use a lumped-system (population-level) modeling (LSM) approach that assumes average properties of a population within a control volume. For modern models that formulate phytoplankton growth as a nonlinear function of the internal nutrient (e.g., Droop kinetics), this averaging assumption can introduce a significant error. Agent-based (individual-based) modeling (ABM) is an alternative approach that does not make the assumption of average properties.

This educational applet contrasts the agent-based and lumped-system modeling approaches. ABM and LSM models are constructed based on identical underlying sub-models of nutrient uptake (Michaelis-Menten) and growth (Droop). The models are applied to a scenario consisting of a point source nutrient discharge into a waterbody with advection and dispersion.

The details of the model can be found in the paper by Hellweger and Kianirad (2007a). More information on ABM vs. LSM comparison can be found in that and the other papers listed below.


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 "Phytoplankton" plot to watch the ABM and LSM populations develop over time.
5. View the Mean and Max. phytoplankton concentrations, and their differences.

Parameters:

VmaxPO4 (molP/molC/d): Max. PO4 uptake rate
KmPO4 (umol/L): PO4 uptake half-saturation constant
uPmax (1/d): Max. photosynthesis rate
q0P (mmolP/molC): P subsitence quota
uR (1/d): Respiration rate
uD (1/d): Death rate
vs (m/d): settling velocity
dt (d): Integration step size
xypatch (m): Grid cell x and y dimension
z (m): Water column depth
vx (m/s): Velocity in x-direction
Eu (m2/s): Dispersion coefficient
WPO4 (mol/s): PO4 loading rate
m0 (pmolC/cell): Min. cell size
rcv (): Randomization coefficent
ChlaC (ugChla/mgC): Chlorophyll a to carbon ratio

Notes:

Several parameters can be changed through the user interface, including VmaxPO4, KmPO4, uPmax, q0P, uR, uD and vs. Other are, for simplicity, are not included in the user interface, but hardwired into the code. This includes dt (0.005 d), xypatch (5,000 m), z (10 m), vx (0.1 m/s), Eu (50 m2/s), WPO4 (0.1 mol/s), m0 (1.0 pmolC/cell), rcv (0.1), ChlaC (30 ugChla/mgC).


THINGS TO NOTICE

Even though the ABM and LSM models are constructed based on identical underlying sub-models, they produce different results. Why is that? Read some of the papers in the reference list to learn more about this.


THINGS TO TRY

Try adjusting the parameters under various settings. How sensitive is the difference between the ABM and LSM models to the particular parameters?


CREDITS AND REFERENCES

Hellweger References:

Hellweger, F. L., Kianirad, E. 2007a. Accounting for Intra-Population Variability in Biogeochemical Models using Agent-Based Methods. Environmental Science & Technology, 41(8):2855-2860.

Hellweger, F. L., Kianirad, E. 2007b. Individual-Based Modeling of Phytoplankton: Evaluating Approaches for Applying the Cell Quota Model. Journal of Theoretical Biology, 249:554-565.

Hellweger, F. L. 2007. Is it time to abandon the chemistry approach to biogeochemistry? (Agent-based water quality modeling). In: Proceedings of WEFTEC 07, Water Environment Federation (WEF), Alexandria, VA., pp. 5646-5665.

Hellweger, F. L. 2008. Spatially Explicit Individual-Based Modeling Using a Fixed Super-Individual Density. Computers & Geosciences, 34(2):144-152.

Hellweger, F. L., Bucci, V. 2009. A bunch of tiny individuals Individual-based modeling for microbes (review paper). Ecological Modelling, 220(1):8-22.

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

; ABMvsLSM.nlogo
; ferdi hellweger
; ferdi@coe.neu.edu
; 12/30/2009

globals 
  [ 
  dt time dtp ntp 
  skipABM 
  xypatch z acpatch vpatch 
  vx vy Qx Qy Eu Ep
  pXa pXb pXt pYa pYb pYt idx idy r1 r2
  xcor2 ycor2
  isx isy WPO4
  vS2
  linear?
  VmaxPO42 KmPO42 VPO4t
  uPmax2 q0P2 uR2 uD2 nq
  m0 SR nb mt sf rcv 
  ChlaC 
  meanLSM meanABM meandiff
  maxLSM maxABM maxdiff
  G1 Cmin Cmax ccolmin ccolmax ccolslope
  LSM_PCt LSM_PPt LSM_PO4t ABM_PO4t
  PC0 PP0 PO40
  PC0s PP0s PO40s
  LSM_RPCt LSM_RPPt LSM_RPO4t ABM_RPO4t
  PO4t
  plotx ploty
  ]

breed 
  [ 
  bguys bguy
  ]

patches-own 
  [ 
  itype 
  LSM_PC LSM_PP LSM_PO4 LSM_TP
  ABM_PC ABM_PP ABM_PO4 ABM_TP
  LSM_RPC LSM_RPP LSM_RPO4 ABM_RPO4
  LSM_VPO4 LSM_qP LSM_uP
  ]

bguys-own 
  [ 
  m VPO4 uP qP RqP
  ]

to setup

  clear-all

  set dt 0.005
  set time 0
  set dtp 1
  set ntp 0 
  set skipABM 0

  set xypatch 5000 
  set z 10
  set acpatch xypatch * z
  set vpatch ( xypatch ^ 2 ) * z

  set vx 0.1 * 86400
  set vy 0
  set Qx vx * acpatch
  set Qy vy * acpatch
  
  set Eu ( 50 * 86400 )
  set Ep ( Eu * acpatch / xypatch )

  set isx -13
  set isy 0
  set WPO4 0.1 * 86400 * ( max-pycor - min-pycor + 1) * xypatch / 1000

  set linear? false
  set m0 1.0e-12
  set rcv 0.1
  if linear? [set KmPO42 KmPO42 / 100]
  set nq 1
  do-input

  set PC0 0
  set PP0 0
  set PO40 0
  set PC0s 1
  set PP0s PC0s * q0P2
  set PO40s KmPO42 * 0.01
  set nb 300

  set SR PC0s * vpatch / (nb * m0 * 1.5)

  set ChlaC 30
  set ccolmax 69
  set ccolmin 61

  ;random-seed 1
  
  ask patches 
    [
    set itype 0
    set LSM_PC PC0
    set LSM_PP PP0
    set LSM_PO4 PO40
    set ABM_PO4 PO40
    ] 
  ask patch isx isy 
    [
    set itype 1
    set LSM_PC PC0s
    set LSM_PP PP0s
    set LSM_PO4 PO40s 
    set ABM_PO4 PO40s 
    set plabel "S "
    set plabel-color black
    ]

  set-default-shape bguys "circle"
  create-bguys nb
    [
    set color 15
    set size 0.15 
    ask one-of patches with [itype = 1] 
      [
      set xcor2 pxcor
      set ycor2 pycor
      ]
    setxy xcor2 ycor2
    set m m0 * 1.5
    set qP PP0s / PC0s
    ]

  do-output
  do-plot

end

to go

  do-input

  if skipABM = 0 
    [
    if not any? turtles [ stop ]
    ]

  if skipABM = 0 [

  ask bguys 
    [

    set idx 0
    set idy 0
    
    set pXb ( Qx / vpatch * dt )
    set r1 random-float 1.0
    if ( r1 < pXb ) [
      set idx 1
      ]

    set pXa ( Ep / vpatch * dt )
    set pXb ( Ep / vpatch * dt )
    set pXt pXa + pXb
    set r1 random-float 1.0
    if ( r1 < pXt ) [
      set r2 random-float 1.0
      ifelse (r2 < pXa / pXt)
        [set idx ( idx - 1 )]
        [set idx ( idx + 1 )]
    ]
    set pYa ( Ep / vpatch * dt )
    set pYb ( Ep / vpatch * dt )
    set pYt pYa + pYb
    set r1 random-float 1.0
    if ( r1 < pYt ) [
      set r2 random-float 1.0
      ifelse (r2 < pYa / pYt)
        [set idy ( idy - 1 )]
        [set idy ( idx + 1 )]
    ]

    set xcor2 ( [pxcor] of patch-here + idx )
    set ycor2 ( [pycor] of patch-here + idy )

    ifelse xcor2 > (min-pxcor - 0.5) and xcor2 < (max-pxcor + 0.5)
    and ycor2 > (min-pycor - 0.5) and ycor2 < (max-pycor + 0.5)
      [ 
      setxy xcor2 ycor2  
      setxy (xcor - 0.5 + random-float 1.0) (ycor - 0.5 + random-float 1.0)
      ]
      [die]
    ]

  ask bguys 
    [
    set PO4t [ABM_PO4] of patch-here 
    set uP 0
    ifelse linear? 
      [
      set uP ( uPmax2 * PO4t / KmPO42 )
      set VPO4 ( (uP - uR2 ) * q0P2 )
      ]
      [
      if qP > q0P2
        [
        set uP ( uPmax2 * ( (qP - q0P2) ^ nq / (q0P2 ^ nq + (qP - q0P2) ^ nq ) ) )
        ]
      set VPO4 ( VmaxPO42 * ( PO4t / ( KmPO42 + PO4t ) ) )
      ]
    set VPO4t VPO4
    set mt m
    ask patch-here [
      set ABM_PO4 ( ABM_PO4 - VPO4t * mt * SR / vpatch * dt )
      ]
    set RqP ( VPO4 - ( uP - uR2 ) * qP)
    set qP ( qP + RqP * dt )
    set m ( m + (uP - uR2) * m * dt )
    if  m > 2 * m0 
      [
      ;type "Division: " type uP type ", " type m type ", " print m0  
      set mt m
      set sf random-normal 0.5 rcv
      if sf <= 0 [set sf 0.5]
      set m ( mt * sf)
      hatch 1 [ 
        set m ( mt * (1 - sf ) )
        ]
      ]
    ]

  ask bguys 
    [
    if random-float 1.0 < ( ( uD2 + vs2 / z) * dt ) [ die ]
    ]

  ] ;skipABM

  ask patches 
    [ 
    set LSM_RPC 0
    set LSM_RPP 0
    set LSM_RPO4 0
    set ABM_RPO4 0    
    ]

  ask patches 
    [ 
    set LSM_PCt [LSM_PC] of self
    set LSM_PPt [LSM_PP] of self
    set LSM_PO4t [LSM_PO4] of self
    set ABM_PO4t [ABM_PO4] of self
    ask neighbors4 
      [
      set LSM_RPC ( LSM_RPC + ( Ep / vpatch * LSM_PCt ) )
      set LSM_RPP ( LSM_RPP + ( Ep / vpatch * LSM_PPt ) )
      set LSM_RPO4 ( LSM_RPO4 + ( Ep / vpatch * LSM_PO4t ) )
      set ABM_RPO4 ( ABM_RPO4 + ( Ep / vpatch * ABM_PO4t ) )
      ]
    set LSM_RPC ( LSM_RPC - ( 4 * Ep / vpatch * LSM_PCt ) )
    set LSM_RPP ( LSM_RPP - ( 4 * Ep / vpatch * LSM_PPt ) )
    set LSM_RPO4 ( LSM_RPO4 - ( 4 * Ep / vpatch * LSM_PO4t ) )
    set ABM_RPO4 ( ABM_RPO4 - ( 4 * Ep / vpatch * ABM_PO4t ) )
    ]

  ask patches 
    [ 
    set LSM_PCt [LSM_PC] of self
    set LSM_PPt [LSM_PP] of self
    set LSM_PO4t [LSM_PO4] of self
    set ABM_PO4t [ABM_PO4] of self
    ask patches at-points [[1 0]] 
      [
      set LSM_RPC ( LSM_RPC + ( Qx / vpatch * LSM_PCt ) )
      set LSM_RPP ( LSM_RPP + ( Qx / vpatch * LSM_PPt ) )
      set LSM_RPO4 ( LSM_RPO4 + ( Qx / vpatch * LSM_PO4t ) )
      set ABM_RPO4 ( ABM_RPO4 + ( Qx / vpatch * ABM_PO4t ) )
      ]
    set LSM_RPC ( LSM_RPC - ( Qx / vpatch * LSM_PCt ) )
    set LSM_RPP ( LSM_RPP - ( Qx / vpatch * LSM_PPt ) )
    set LSM_RPO4 ( LSM_RPO4 - ( Qx / vpatch * LSM_PO4t ) )
    set ABM_RPO4 ( ABM_RPO4 - ( Qx / vpatch * ABM_PO4t ) )
    ]

  ask patches with [itype = 1] 
    [
    set LSM_RPO4 ( LSM_RPO4 + WPO4 / vpatch )
    set ABM_RPO4 ( ABM_RPO4 + WPO4 / vpatch )
    ]

  ask patches 
    [
    set LSM_VPO4 ( VmaxPO42 * ( LSM_PO4 / ( KmPO42 + LSM_PO4 ) ) )
    set LSM_qP 0 
    if LSM_PC > 0 
      [ 
      set LSM_qP ( LSM_PP / LSM_PC ) 
      set LSM_uP 0
      ifelse linear?
        [
        set LSM_uP ( uPmax2 * LSM_PO4 / KmPO42 )
        set LSM_VPO4 ( (LSM_uP - uR2) * q0P2 )
        ]
        [
        if LSM_qP > q0P2
          [
          set LSM_uP ( uPmax2 * ( (LSM_qP - q0P2) ^ nq / (q0P2 ^ nq + (LSM_qP - q0P2) ^ nq ) ) )
          ] 
        ]
      set LSM_RPC ( LSM_RPC + (LSM_uP - uR2 - uD2 - vs2 / z) * LSM_PC )
      set LSM_RPP ( LSM_RPP + LSM_VPO4 * LSM_PC + (- uD2 - vs2 / z) * LSM_PP )
      set LSM_RPO4 ( ( LSM_RPO4 - LSM_VPO4 * LSM_PC ) )
      ]
    ]

  ask patches 
    [
    set LSM_PC ( LSM_PC + LSM_RPC * dt) 
    set LSM_PP ( LSM_PP + LSM_RPP * dt ) 
    set LSM_PO4 ( LSM_PO4 + LSM_RPO4 * dt ) 
    set ABM_PO4 ( ABM_PO4 + ABM_RPO4 * dt ) 
    ]

  if min [LSM_PC] of patches < 0 [show "Error: Negative LSM_PC."]
  if min [LSM_PP] of patches < 0 [show "Error: Negative LSM_PP."]
  if min [LSM_PO4] of patches < 0 [show "Error: Negative LSM_PO4."]

  if ntp >= time 
    [
    set ntp time + dtp
    do-output
    do-plot
    ]

  set time ( time + dt )

end

to do-input

  set VmaxPO42 VmaxPO4 ; 0.5
  set KmPO42 KmPO4 * 1e-3 ; 2.5 * 1e-3
  set uPmax2 uPmax ; 2.0
  set q0P2 q0P * 1e-3 ; 1e-3
  set uR2 uR ; 0.1
  set uD2 uD ; 0.15
  set vs2 vs ;0.5
 
end

to do-output

  ask patches 
    [
    set LSM_TP ( LSM_PP + LSM_PO4 )
    ] 

  set Cmin min [LSM_PC] of patches
  set Cmax max [LSM_PC] of patches
  set ccolslope 0
  if not (Cmax = 0) 
    [
    set ccolslope (ccolmax - ccolmin) / ( Cmax - Cmin) 
    ]
  ask patches 
    [
    set G1 ( Ccolmax -  ( LSM_PC - Cmin ) * ccolslope )
    set pcolor G1
    ]

  ask patches 
    [
    set ABM_PC 0
    set ABM_PP 0
    if any? bguys-here 
      [
      set ABM_PC ( sum [m] of bguys-here ) * SR / vpatch
      set ABM_PP ( ABM_PC * mean [qP] of bguys-here )
      ]
    set ABM_TP ( ABM_PP + ABM_PO4 ) 
    ]

  set meanLSM mean [LSM_PC] of patches * 12 * ChlaC
  set meanABM mean [ABM_PC] of patches * 12 * ChlaC
  ifelse ( ( meanLSM + meanABM ) > 0 )
    [set meandiff abs (meanLSM - meanABM) / ((meanLSM + meanABM) / 2 ) * 100]
    [set meandiff 0]

  set maxLSM max [LSM_PC] of patches * 12 * ChlaC
  set maxABM max [ABM_PC] of patches * 12 * ChlaC
  ifelse ( ( meanLSM + meanABM ) > 0 )
    [set maxdiff abs (maxLSM - maxABM) / ((maxLSM + maxABM) / 2 ) * 100]
    [set maxdiff 0]

end


to do-plot

  set-current-plot "Nutrient"
  set-current-plot-pen "LSMTP"
  plotxy time mean [LSM_TP] of patches * 1e3
  set-current-plot-pen "LSMPP"
  plotxy time mean [LSM_PP] of patches * 1e3
  set-current-plot-pen "LSMPO4"
  plotxy time mean [LSM_PO4] of patches * 1e3
  set-current-plot-pen "ABMTP"
  plotxy time mean [ABM_TP] of patches * 1e3
  set-current-plot-pen "ABMPP"
  plotxy time mean [ABM_PP] of patches * 1e3
  set-current-plot-pen "ABMPO4"
  plotxy time mean [ABM_PO4] of patches * 1e3

  set-current-plot "Phytoplankton"
  set-current-plot-pen "LSM"
  plotxy time meanLSM
  set-current-plot-pen "ABM"
  plotxy time meanABM

end