QuotaHeterogeneity: Modeling Intracellular P Heterogeneity in Cultured Phytoplankton

(c) Neil Fredrick, Ferdi Hellweger, Northeastern University, Boston.

NOTE: This 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.


powered by NetLogo

view/download model file: QuotaHeterogeneity.nlogo

WHAT IS IT?

This model was created to explore mechanisms of heterogeneity for internal P content (mol P/cell) of phytoplankton in laboratory culture. There can be significant intraspecific individual-level heterogeneity in the internal P content of phytoplankton, which can affect the population-level growth rate. Several mechanisms can create this heterogeneity, including biological variability in states and behavior. Here we use modeling to explore the contribution of various mechanisms to the heterogeneity in phytoplankton grown in a laboratory culture. An agent-based model simulates individual cells and their nutrient content. Heterogeneity is introduced by randomizing states (e.g., P content) and parameters (e.g., maximum uptake rate) of daughter cells at division. The model was calibrated to observations of internal P content of individual cells of Thalassiosira pseudonana obtained using synchrotron X-ray fluorescence (SXRF). Then, a number of simulations with individual mechanisms of heterogeneity turned off were performed. Comparison of the coefficient of variation (CV) of these and the baseline simulation (all mechanisms turned on) provides an estimate of the relative contribution of these mechanisms. The results show that the mechanism with the largest contribution is variability in maximum internal cell quota parameter, which when removed results in a CV of 0.21 compared to a CV of 0.38 with all mechanisms turned on. This suggests that nutrient/element storage capabilities/mechanisms are important determinants of intrapopulation variability.

HOW IT WORKS

The ABM simulates one type of phytoplankton (Thalassiosira pseudonana) in culture. The reactor is assumed to be completely mixed (uniform concentration of extracellular nutrient), so the agents do not see any spatial variation in the environment. Cell growth is a function of photosynthesis and respiration. Photosynthesis depends on the internal P content and light, with no effect from temperature and other nutrients. Phytoplankton take up and excrete P. The model includes mortality, and dead cells are tracked. They excrete nutrient but do not perform uptake, photosynthesis, respiration or growth. The details of the model can be found in the paper Fredrick et al (2012).

HOW TO USE IT

Click the SETUP button to create the starting population of phytoplankton cells. Click the GO bottom to start the simulation. The cells change from black to green depending on their internal P content.

There are different mechanisms that create heterogeneity, these mechanisms are controlled by 12 switches that each turn a mechanism on or off.

Mechanisms
Upmax(1/day): Maximum photosynthesis rate
q0(mol P/mol C): P subsistence quota
uR(1/day): Respiration rate
m0(mol C/cell): Minimum cell size
uDbase(1/day): Base death rate
qD(mol P/mol C): P death quota
VAmax(mol P/um2/day): Maximum specific uptake rate
Km(mol P/m3): Half-saturation constant
qmax(mol P/mol C): P maximum quota
kw(1/day): Excretion constant
Sf-qm: Split fraction of internal P
Sf-m: Split fraction of internal C

There is a cumulative density function (CDF) plot for internal P content showing the model output and the laboratory data. Also three plots that demonstrate the differences between the ABM and an equivalent PLM model. They show the differences in extracellular P, cell quota, and the number of cells (living dead and total).

There are also monitors for the time in days, Root mean squared error (RMSE), and model CV.

THINGS TO NOTICE

The all cells in the ABM start with the same mass and quota and it takes some time (about 1.5 days) for the cells to desynchronize. Despite have the same average parameters values the populations in the ABM and PLM the average growth rates of the populations are not the same.

THINGS TO TRY

Try turning off/on various mechanisms (some mechanism have a larger effects on the heterogeneity of the population than others). Also try starting the model with no mechanisms turned on then turn on mechanisms individually to see their effect on heterogeneity.

EXTENDING THE MODEL

This version of the model is for distribution so there are some modifications from the original. Modifications include: The models is setup in chemostat mode with a flow rate equal to the average growth rate. (2) The light is continuous (3)The range for the total number of agents are reduced from 1100-2200 to 550-1100 (4) random seed was set to 1, (5) the mechanisms are turned on/off with switches. If the modifications are removed the model can be run as a batch culture or a periodically diluted batch culture.

There is also code for an optimization routine that finds the best CV’s for the mechanisms that are switched on/off in this version. The optimization routine and its associated variables have been commented out with “;;OPT::” to make the code simpler and therefore easier to understand.

RELATED MODELS

ABMvsLSM: Agent-Based Model vs. Lumped-System Model (v1.0)

CREDITS AND REFERENCES

Fredrick, N. D., Berges, J. A., Twining B. S., Nuñez-Milland, D., Hellweger, F. L. (2012). Exploring Mechanisms of P Content Heterogeneity in Cultured Phytoplankton using Agent Based Modeling. In Press

CODE

;Notes:
;1.  This model corresponds to the M1 simulation described in the accompanying
; paper.  However, this distribution version of the model was modified from the one
;used for the simulations presented in the accompanying paper to factilitate
;learning and adoption by new users.  Modifications include: (a) The models is setup
;in chemostat mode with a flow rate equal to the average growth rate (b) light is 
;continuous (c) the range for the total number of agents are reduced from 1100-2200
;to 550-1100 (d) the random seed was set to 1, (e) the mechanisms are turned on/off
;with switches.  The optimization routine and its associated variables have been 
;commented out with ";;OPT::" to make the code simpler and therefore easier to 
;understand.
; 2.  The model includes an optimization routine (located at the end of the code).
; The routine changes the CV's for some of the model parameters (uPmax, q0, uR, m0,
; uDbase, qD, VAmax, Km, qmax, kw, Sf-qm, and Sf-m).  To change the parametes CV's
; every parameter has a base CV (i.e. uPmax-CV-base) and a modifying variable 
; (i.e. var-1) which is what is varied during optimization.  The base is multiplied
; by the modifying variable to get the CV.  The optimization routine changes the
; modifying variables.
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;Name: iAlgae
;Version: n.nl.1.0
;Date: March 21, 2012
;Authors: Neil Fredrick(1), John Berges(2) and Ferdi Hellweger(1)
;     (1) Department of Civil and Environmental Engineering
;         Northeastern University
;     (2) Department of Biological Sciences and School of Freshwater Science
;         University of Wisconsin-Milwaukee
;Desciption: ABM model of phytoplankton in a laboratory reactor
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Global Variables;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

globals [
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Phytoplankton;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Parameters;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 uPmax-Mean                       ;max photosynthesis rate (1/d)-Mean
 q0-Mean                          ;subsistence quota (molP/molC)-Mean
 uR-Mean                          ;respiration rate (1/d)-Mean
 m0-Mean                          ;minimum cell size (molC/cell)-Mean
 uDbase-Mean                       ;maximum death rate (1/d)-Mean
 qD-Mean                          ;death quota (molP/molC)-Mean
 VAmax-Mean                       ;maximum uptake rate (molP/um2/d)-Mean
 Km-Mean                          ;half-saturation constant (molP/m3)-Mean
 qmax-Mean                        ;maximum quota (molP/molC)-Mean
 kw-Mean                          ;excretion constant (1/d)-Mean
 Sf-qm-Mean                       ;split-fraction for quota-Mean
 Sf-m-Mean                        ;split-fraction for size-Mean
 
 uPmax-CV                         ;CV for uPmax
 q0-CV                            ;CV for q0
 uR-CV                            ;CV for uR
 m0-CV                            ;CV for m0
 uDbase-CV                         ;CV for uDbase
 qD-CV                            ;CV for qD
 VAmax-CV                         ;CV for VAmax
 Km-CV                            ;CV for Km
 qmax-CV                          ;CV for qmax
 kw-CV                            ;CV for kw
 Sf-qm-CV                         ;CV for Sf-qm
 Sf-m-CV                          ;CV for Sf-m
 
 p                                ;cell density (molC/um3)
 d                                ;cell diameter (um)
 C0                               ;initial cell concentration(cells/ml)
 
 Data-x                           ;list of data (quotas)
 Data-y                           ;list of data (probabilities)
 m-plot-scale                     ;size plotting scale
 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Calculated Values;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 
 Sf-qm                            ;split-fraction for quota
 Sf-m                             ;split-fraction for size

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Run Control;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Parameters;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 dt                               ;time step (d)
 dtp                              ;time step-print (d)
 tend                             ;end of simulation time
 ;randomseed                       ;random seed
 file-prefix                      ;prefix output file
 file-suffix                      ;suffix output file
 
 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;State Variables;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 time                             ;time (d)
 file-number                      ;output file number
 filename                         ;output file name 1
 filename-2                       ;output file name 2
 

 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Agent Accounting;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Parameters;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 n-min-alive                      ;min number of "Super Individuals" alive
 n-max-alive                      ;max number of "Super Individuals" alive
 n-min-dead                       ;min number of "Super Individuals" dead
 n-max-dead                       ;max number of "Super Individuals" dead
 n-start                          ;starting number "Super Individuals" alive
 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Calculated Values;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 Sr-min                           ;Super Individual number of a cell being combined
 m-min                            ;size of a cell being combined
 age-min                          ;age of a cell being combined
 q-min                            ;quota of a cell being combined
 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Environment;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Parameters;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 Vol                              ;volume (m3)
 S0                               ;initial concentration of PO4 in reactor (molP/m3)
 S-in                             ;concentration of PO4 inflow (molP/m3)
 Flow                             ;flow rate (m3/d)
 dt-d                             ;dilution time step (d)
 dtL                              ;time of light (d)
 dtD                              ;time of Dark (d)
 f-d                              ;dilution ratio

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;State Variables;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 S                                ;concentration of PO4 in reactor (mol/m3)
 LI                               ;light limitation term (1 when light, 0 when dark)
 tlnext                           ;next time to change LightDark value (d)
 tdnext                           ;next dilution time (d)
 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Calculated Values;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

 sumS                             ;total P in the reactor (intra and extra cellular)
 sumq                             ;P in cells in reactor
 ave-q-dead                       ;average quota of dead cells
 weighted-average-q               ;average that takes Sr into account
 weighted-average-m
 weighted-average-qm              ;average of q*m that takes Sr into account
 weighted-std-dev-q               ;std-dev that takes Sr into account
 weighted-std-dev-m
 weighted-CV                      ;CV that takes Sr into account
 weighted-VAR                     ;variance of q*m that takes Sr into account
 weighted-Var-m
 RMSE                             ;Root Mean Squared Error 

;;OPT::;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;OPT::;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Optimization;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;OPT::;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;OPT:: 
;;OPT::;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Parameters;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;OPT:: t-end-opt                        ;time for run in optimization
;;OPT:: uPmax-CV-base                    ;CV for uPmax-base for optimization
;;OPT:: q0-CV-base                       ;CV for q0-base for optimization
;;OPT:: uR-CV-base                       ;CV for uR-base for optimization
;;OPT:: m0-CV-base                       ;CV for m0-base for optimization
;;OPT:: uDbase-CV-base                   ;CV for uDbase-base for optimization
;;OPT:: qD-CV-base                       ;CV for qD-base for optimization
;;OPT:: VAmax-CV-base                    ;CV for VAmax-base for optimization
;;OPT:: Km-CV-base                       ;CV for Km-base for optimization
;;OPT:: qmax-CV-base                     ;CV for qmax-base for optimization
;;OPT:: Kw-CV-base                       ;CV for Kw-base for optimization
;;OPT:: Sf-m-CV-base                     ;CV for Sf-m-base for optimization
;;OPT:: Sf-qm-CV-base                    ;CV for Sf-qm-base for optimization
;;OPT:: var-1                            ;number on variable list for a variable
;;OPT:: var-2                            ;number on variable list for a variable
;;OPT:: var-3                            ;number on variable list for a variable
;;OPT:: var-4                            ;number on variable list for a variable
;;OPT:: var-5                            ;number on variable list for a variable
;;OPT:: var-6                            ;number on variable list for a variable
;;OPT:: var-7                            ;number on variable list for a variable
;;OPT:: var-8                            ;number on variable list for a variable
;;OPT:: var-9                            ;number on variable list for a variable
;;OPT:: var-10                           ;number on variable list for a variable
;;OPT:: var-11                           ;number on variable list for a variable
;;OPT:: var-12                           ;number on variable list for a variable
;;OPT:: shuffle-order                    ;list of numbers to randomly assigned to a variable
;;OPT:: opt-step-size                    ;size of percent change in optimization
;;OPT:: opt-runs-step                    ;number of runs for the error to be averaged over
;;OPT:: 
;;OPT::;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;State Variables;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;OPT:: opt-fact                         ;list of CV for variables
;;OPT:: var-new                          ;new variable value
;;OPT:: var-current                      ;starting variable value
;;OPT:: RMSE-min                         ;starting root mean squared error
;;OPT:: opt-run-number                   ;number of runs
;;OPT:: i                                ;looping variable
;;OPT:: j                                ;looping variable
;;OPT:: k                                ;looping variable
;;OPT:: opt-step-total                   ;total error for current step

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;PLM Variables;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Parameters;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 PLM-uPmax                        ;max photosynthesis rate (1/d)-PLM
 PLM-q0                           ;subsistence quota (molP/molC)-PLM
 PLM-ur                           ;repiration rate (1/d)-PLM
 PLM-uDbase                       ;maximum death rate (1/d)-PLM
 PLM-qD                           ;death quota (molP/molC)-PLM
 PLM-Vmax                         ;maximum uptake rate (molP/um2/d)-PLM
 PLM-Km                           ;half-saturation constant (molP/m3)-PLM
 PLM-qmax                         ;maximum quota (molP/molC)-PLM
 PLM-kw                           ;excretion constant (1/d)-PLM
 PLM-mean-cell-C                  ;average phytoplankton C-PLM
 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;State Variables;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 PLM-S                            ;concentration of PO4 in reactor (mol/m3)-PLM
 PLM-ug                           ;growth rate-PLM
 PLM-V                            ;uptake rate-PLM
 PLM-W                            ;excretion rate-PLM
 PLM-q                            ;quota-PLM
 PLM-uD                           ;mortality rate(1/d)-PLM
 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Calculated Values;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 PLM-cell-C                       ;phytoplankton C-PLM
 PLM-dead-cell-C                  ;dead phytoplankton C-PLM
]

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Agent Variables;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
                                  
turtles-own [
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Phytoplankton;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Parameters;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 uPmax                            ;maximum photosynthesis rate (1/d)
 q0                               ;subsistence cell quota (molP/molC)
 uR                               ;respiration rate (1/d)
 m0                               ;minimum cell size (molC/cell)
 uDbase                           ;max death rate (1/d)
 qD                               ;death quota (molP/molC)
 VAmax                            ;maximum uptake rate area (molP/um2/day)
 Km                               ;half-saturation constant (molP/m3)
 qmax                             ;maximum cell quota (molP/molC)
 kw                               ;excretion constant (1/d)
 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;State Variables;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 m                                ;cell size (molC/cell)
 q                                ;cell quota (molP/molC) 
 Sr                               ;"Super Individual" representative number
 age                              ;cell age (d)
 alive                            ;if the cell is alive or dead (1 or 0)
 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Calculated Values;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
 uP                               ;photosynthesis rate (1/d)
 uG                               ;growth rate (1/d)
 uD                               ;death rate (1/d)
 q-m                              ;cell quota (molP/cell) 
 V                                ;uptake rate (molP/molC/d)
 W                                ;excretion rate (molP/molC/d)
 q-Sr                             ;quota * Sr
 q-m-Sr                           ;quota * size * Sr
 q-w                              ;weighted quota
 cell-Vol                         ;volume of cell (um3)
 A                                ;cell surface area (um2)
 h                                ;cell height (um)
 VWm-Sr                           ;net uptake of super individual (molP/d)
 Pm                               ;probability a cell will die this time step
 CDF-probability                  ;agents probability in the CDF plot
]

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Setup;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

To setup
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Run Control;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Initialize Netlogo;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  clear-all
  reset-ticks
  random-seed 1

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Parameters;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  set dt .0005
  set dtp 1
  set tend 7.5
  set m-plot-scale 10 ^ 12.5
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;State Variables;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  set time 0
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Agent Accounting;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Parameters;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  set n-min-alive 500
  set n-max-alive 1000
  set n-min-dead 50
  set n-max-dead 100
  set n-start (n-min-alive + n-max-alive) / 2
  
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Environment;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Parameters;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  set Vol 0.001
  set S0 0.0015
  set S-in S0
  set Flow .0013
  set dt-d 999999
  set dtL 99999999
  set dtD 10 / 24
  set f-d .99
 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;State Variables;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  set S S0
  set LI 1
  set tdnext dt-d
  set tlnext dtL
  ask patches [set pcolor white]

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Phytoplankton;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Parameters;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;  
  set uPmax-Mean 1.6
  set q0-Mean 1.52e-3
  set uR-Mean 0.1
  set m0-Mean 1.8e-13
  set uDbase-Mean 0.20
  set qD-Mean 1.8e-3
  set VAmax-Mean 2.64e-15
  set Km-Mean 6e-4
  set qmax-Mean 1.5e-2
  set kw-Mean 0.01
  set Sf-m-Mean .5
  set Sf-qm-Mean .5
;;OPT::  
;;OPT::  set uPmax-CV-base 0.1
;;OPT::  set q0-CV-base 0.1
;;OPT::  set uR-CV-base 0.1
;;OPT::  set m0-CV-base 0.1
;;OPT::  set uDbase-CV-base 0.1
;;OPT::  set qD-CV-base 0.1
;;OPT::  set VAmax-CV-base 0.1
;;OPT::  set Km-CV-base 0.1
;;OPT::  set qmax-CV-base 0.1
;;OPT::  set Kw-CV-base 0.1
;;OPT::  set Sf-qm-CV-base 0.1
;;OPT::  set Sf-m-CV-base 0.1
  
  set p 2.46 * 10 ^ -14
  set d 4.13
  set C0 50000 * 10 ^ 6 ;converts cells/ml to cells/m3
  
  set Data-x [1.20891562005429 1.41288793958447 1.47342101642315 1.79431241925361 1.92575030758369 2.0737622827365 2.10831985694315 2.55018571053791 2.6679602295747 2.79128032987699 2.91349282594903 3.23057444320554 3.27656215610311 3.48929825767127 3.52087216841275 3.60035397933587 3.74111568703057 3.872373154644 4.16968290464708 4.23505571489438 4.24918639945223 4.25937147264245 4.29028990337407 4.34855060063512 4.39024895051015 4.43257453847461 4.43358057405921 4.46814721027696 4.47756202511576 4.48673759554141 4.63400629820206 4.88901052527293 4.91072100177276 5.17962994276562 5.19488507136414 5.27886151729539 5.45109544169034 5.62848815448372 6.19139051749515 6.21812337545462 6.70736617935475]
  set Data-y [0.0238095238095238 0.0476190476190476 0.0714285714285714 0.0952380952380952 0.119047619047619 0.142857142857143 0.166666666666667 0.19047619047619 0.214285714285714 0.238095238095238 0.261904761904762 0.285714285714286 0.30952380952381 0.333333333333333 0.357142857142857 0.380952380952381 0.404761904761905 0.428571428571429 0.452380952380952 0.476190476190476 0.5 0.523809523809524 0.547619047619048 0.571428571428571 0.595238095238095 0.619047619047619 0.642857142857143 0.666666666666667 0.69047619047619 0.714285714285714 0.738095238095238 0.761904761904762 0.785714285714286 0.80952380952381 0.833333333333333 0.857142857142857 0.880952380952381 0.904761904761905 0.928571428571429 0.952380952380952 0.976190476190476]
  
;;OPT::;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;OPT::;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;optimization;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;OPT::;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;  
;;OPT::  
;;OPT::  set shuffle-order [0 1 2 3 4 5 6 7 8 9 10 11]
;;OPT::  set opt-fact [0 1.5 2 0 1 1.5 1 .5 6 0 4 .5]
;;OPT::  set var-1 item 0 shuffle-order
;;OPT::  set var-2 item 1 shuffle-order
;;OPT::  set var-3 item 2 shuffle-order
;;OPT::  set var-4 item 3 shuffle-order
;;OPT::  set var-5 item 4 shuffle-order
;;OPT::  set var-6 item 5 shuffle-order
;;OPT::  set var-7 item 6 shuffle-order
;;OPT::  set var-8 item 7 shuffle-order
;;OPT::  set var-9 item 8 shuffle-order
;;OPT::  set var-10 item 9 shuffle-order
;;OPT::  set var-11 item 10 shuffle-order
;;OPT::  set var-12 item 11 shuffle-order
;;OPT::
;;OPT::  set uPmax-CV (uPmax-CV-base * item var-1 opt-fact)
;;OPT::  set q0-CV (q0-CV-base * item var-2 opt-fact)
;;OPT::  set uR-CV (uR-CV-base * item var-3 opt-fact)
;;OPT::  set m0-CV (m0-CV-base * item var-4 opt-fact)
;;OPT::  set uDbase-CV (uDbase-CV-base * item var-5 opt-fact) 
;;OPT::  set qD-CV (qD-CV-Base * item var-6 opt-fact)
;;OPT::  set VAmax-CV (VAmax-CV-base * item var-7 opt-fact)
;;OPT::  set Km-CV (Km-CV-base * item var-8 opt-fact)
;;OPT::  set qmax-CV (qmax-CV-base * item var-9 opt-fact)
;;OPT::  set Kw-CV (Kw-CV-base * item var-10 opt-fact)
;;OPT::  set Sf-qm-CV (Sf-qm-CV-base * item var-11 opt-fact)
;;OPT::  set Sf-m-CV (Sf-m-CV-base * item var-12 opt-fact)

ifelse Upmax-Variability
[set uPmax-CV 0]
[set uPmax-CV 0]
ifelse q0-Variability
[set q0-CV 0.15]
[set q0-CV 0]
ifelse uR-Variability
[set uR-CV 0.2]
[set uR-CV 0]
ifelse m0-Variability
[set m0-CV 0]
[set m0-CV 0]
ifelse uDbase-Variability
[set uDbase-CV 0.1]
[set uDbase-CV 0]
ifelse qD-Variability
[set qD-CV 0.15]
[set qD-CV 0]
ifelse VAmax-Variability
[set VAmax-CV 0.1]
[set VAmax-CV 0]
ifelse Km-Variability
[set Km-CV 0.05]
[set Km-CV 0]
ifelse qmax-Variability
[set qmax-CV 0.6]
[set qmax-CV 0]
ifelse kw-Variability
[set kw-CV 0]
[set kw-CV 0]
ifelse Sf-qm-Variability
[set Sf-qm-CV 0.4]
[set Sf-qm-CV 0]
ifelse Sf-m-Variability
[set Sf-m-CV 0.05]
[set Sf-m-CV 0]

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Agent Creation;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 
  crt n-start
    [
    set uPmax kreft uPmax-Mean uPmax-CV 0 99
    set q0 kreft q0-Mean q0-CV 0 99
    set uR kreft uR-Mean uR-CV 0 99
    set m0 kreft m0-Mean m0-CV 0 99
    set uDbase kreft uDbase-Mean uDbase-CV 0 99
    set qD kreft qD-Mean qD-CV 0 99
    set VAmax kreft VAmax-Mean VAmax-CV 0 99
    set Km kreft Km-Mean Km-CV 0 99
    set qmax kreft qmax-Mean qmax-CV 0 99
    set kw kreft kw-Mean kw-CV 0 99
    set m m0-Mean
    set q  q0-Mean * 5
    set Sr C0 * vol / n-start
    set age 0
    set alive 1 
    set h m / p * 4 / (pi * d ^ 2)
    set A .5 * pi * d ^ 2 + pi * d * h
    set CDF-probability 0
    set size m * m-plot-scale
    set shape "circle"
    set color green
    set xcor random-xcor
    set ycor random-ycor
    ]
    move-turtles-1
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;PLM setup;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;     
 set PLM-S S0
 set PLM-mean-cell-C m0-Mean;2.49e-13
 set PLM-cell-C C0 * vol * PLM-mean-cell-C
 set PLM-uPmax uPmax-Mean
 set PLM-q0 q0-Mean
 set PLM-ur uR-Mean
 set PLM-uDbase uDbase-Mean
 set PLM-qD qD-Mean
 set PLM-Vmax .388                         ;maximum uptake rate (molP/um2/d)
 set PLM-Km Km-Mean
 set PLM-qmax qmax-Mean
 set PLM-kw kw-Mean
 ;set PLM-V                            ;uptake rate-PLM
 ;PLM-W                            ;excretion rate-PLM
 set PLM-q q0-Mean * 5
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Main Go Loop;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

To go
  tick
  ask turtles with [alive = 1] [uptake-rate]
  ask turtles [excretion-rate]
  ask turtles [mass-balance-cells]
  mass-balance-chemostat
  ask turtles with [alive = 1] [growth-rate]
  ask turtles with [alive = 1] [grow]
  ask turtles with [alive = 1] [update-age]
  ask turtles with [alive = 1] [divide]
  ask turtles with [alive = 1] [mortality]
  if flow > 0 [ask turtles [wash-out]]
  if time > tdnext
  [dilute]
  if count turtles < 1 [stop]
  split-combine
  move-turtles-1
  if time > tlnext [update-light]
  set time time + dt
  update-plot
  update-plot-CDF
  calculated-stats
  Root-Mean-Squared-Error
  PLM
  if time > tend [update-plot-CDF calculated-stats stop]
  check-Variability
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Uptake and Excretion Routines;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
to uptake-rate
   set A .5 * pi * d ^ 2 + pi * d * h
   ifelse (S > 0) and (q < qmax) and (q0 < qmax)
   [set V VAmax * A / m * S / (Km + S) * (qmax - q) / (qmax - q0)]
   [set V 0]
end

to excretion-rate
  ifelse (q > 0)
    [set W kw * q]
    [set W 0]
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Mass Balances;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
to mass-balance-cells
  set VWm-Sr (V - W) * m * Sr
  set q q + (V - W - q * ug) * dt
  if q < 0 [stop]
end

to mass-balance-chemostat
  set S S + (Flow / Vol * S-in - Flow / Vol * S - sum [VWm-Sr] of turtles / Vol) * dt
  ; if s < 0 [set s 0]
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Growth Routines;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
to growth-rate
  ifelse q > q0 
    [set uP uPmax * LI * (1 - q0 / q)]
    [set uP 0]
    set uG uP - uR
end

to grow
  set m m + m * uG * dt
  set h m / p * 4 / (pi * d ^ 2)
  ;set cell-Vol  pi * (d / 2) ^ 2 * h
  set size m * m-plot-scale
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Age Routine;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
to update-age
  set age age + dt
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Cell Division Routine;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
to divide
  if m >= 2 * m0 [
    let Mm m
    let Mqm q * m
    set age 0
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Randomize States;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    set Sf-m kreft Sf-m-Mean Sf-m-CV 0.1 0.9
    set Sf-qm kreft Sf-qm-Mean Sf-qm-CV 0.1 0.9  
    set m Mm * Sf-m
    set q Mqm * Sf-qm / m
    set h m / p * 4 / (pi * d ^ 2)
      
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Randomize Parameters;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    set uPmax kreft uPmax-Mean uPmax-CV 0 99
    set q0 kreft q0-Mean q0-CV 0 99
    set uR kreft uR-Mean uR-CV 0 99
    set m0 kreft m0-Mean m0-CV 0 99
    set uDbase kreft uDbase-Mean uDbase-CV 0 99
    set qD kreft qD-Mean qD-CV 0 99
    set VAmax kreft VAmax-Mean VAmax-CV 0 99
    set Km kreft Km-Mean Km-CV 0 99
    set qmax kreft qmax-Mean qmax-CV 0 99
    set kw kreft kw-Mean kw-CV 0 99

     hatch 1 [
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Randomize States;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
       set m Mm * (1 - Sf-m)
       set age 0
       set q Mqm * (1 - Sf-qm) / m
       set h m / p * 4 / (pi * d ^ 2)
       
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Randomize Parameters;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    set uPmax kreft uPmax-Mean uPmax-CV 0 99
    set q0 kreft q0-Mean q0-CV 0 99
    set uR kreft uR-Mean uR-CV 0 99
    set m0 kreft m0-Mean m0-CV 0 99
    set uDbase kreft uDbase-Mean uDbase-CV 0 99
    set qD kreft qD-Mean qD-CV 0 99
    set VAmax kreft VAmax-Mean VAmax-CV 0 99
    set Km kreft Km-Mean Km-CV 0 99
    set qmax kreft qmax-Mean qmax-CV 0 99
    set kw kreft kw-Mean kw-CV 0 99
       ]
     ]
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Mortality Routines;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
to mortality
  set uD uDbase * qD / q
  set Pm uD * dt
  if Pm > random-float 1
  [
    set alive 0
    set V 0
    set uR 0
    set uG 0
    set shape "square"
  ]
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Washout;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
to wash-out
  set Sr Sr - Flow / Vol * Sr * dt
  if Sr <= 1 [die]
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Dilution Routine;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
to dilute
  set tdnext time + dt-d
  ask turtles [
    set Sr Sr * (1 - f-d)
    if Sr < 1 [die]
    ]
  set PLM-cell-C PLM-cell-C * (1 - f-d)
  set PLM-dead-cell-C PLM-dead-cell-C * (1 - f-d)
  
  set S S0 * f-d + S * (1 - f-d)
  
  set PLM-S S0 * f-d + PLM-S * (1 - f-d)
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Light/Dark Routine;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
to update-light
  if (time >= tlnext) and (LI = 1)
  [
    set LI 0
    set tlnext tlnext + dtD 
  ]
  if (time >= tlnext) and (LI = 0)
  [
    set LI 1
    set tlnext tlnext + dtL
  ]
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Kreft Function;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
to-report kreft [average cv low-limit high-limit]
  ifelse cv = 0
  [
    report average
   ]
  [
  let stdev average * cv
  let temp 0
 while [(temp > (average + 2 * stdev)) or (temp < (average - 2 * stdev)) or (temp <= low-limit) or (temp >= high-limit)]
  [set temp random-normal average stdev]
  report temp
    ]
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Split/Combine Routine;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
to split-combine
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Combine Alive;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
while [count turtles with [alive = 1] > n-max-alive]
   [
    ask min-one-of turtles with [alive = 1] [Sr]
    [
      set Sr-min Sr
      set m-min m
      set age-min age
      set q-min q
      die
      ]
    ask min-one-of turtles with [alive = 1] [Sr]
    [
      set q (q-min * m-min * Sr-min + q * m * Sr) / (Sr-min + Sr)
      set m (m-min * Sr-min + m * Sr) / (Sr-min + Sr)
      set q q / m
      set age (age-min * Sr-min + age * Sr) / (Sr-min + Sr)
      set Sr Sr + Sr-min
     ]
    ]
   
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Split Alive;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
while [count turtles with [alive = 1] < n-min-alive and max [Sr] of turtles with [alive = 1]  > 1]
     [
    ask max-one-of turtles with [alive = 1] [Sr]
    [
      set Sr Sr / 2
      hatch 1
      ]
    ]
     
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Combine Dead;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
if count turtles with [alive = 0] > 0 [
while [count turtles with [alive = 0] > n-max-dead]
   [
    ask min-one-of turtles with [alive = 0] [Sr]
    [
      set Sr-min Sr
      set m-min m
      set age-min age
      set q-min q
      die
      ]
    ask min-one-of turtles with [alive = 0] [Sr]
    [
      set q (q-min * m-min * Sr-min + q * m * Sr) / (Sr-min + Sr)
      set m (m-min * Sr-min + m * Sr) / (Sr-min + Sr)
      set q q / m
      set age (age-min * Sr-min + age * Sr) / (Sr-min + Sr)
      set Sr Sr + Sr-min
     ]
    ]
   ]

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Split Dead;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
if count turtles with [alive = 0] > 0 [       
while [(count turtles with [alive = 0] < n-min-dead) and max [Sr] of turtles with [alive = 0] > 1]
     [
    ask max-one-of turtles with [alive = 0] [Sr]
    [ 
    set Sr Sr / 2
      hatch 1
      ]
     ]
    ]
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Plotting Routines;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
to update-plot
  set-current-plot "Cells"
  set-current-plot-pen "Alive"
  plotxy time (sum [Sr] of turtles with [alive = 1])
  set-current-plot-pen "Dead"
  plotxy time (sum [Sr] of turtles with [alive = 0])
  set-current-plot-pen "Total"
  plotxy time (sum [Sr] of turtles)
  set-current-plot-pen "PLM-Alive"
  plotxy time (PLM-cell-C / PLM-mean-cell-C)
  set-current-plot-pen "PLM-Dead"
  plotxy time (PLM-dead-cell-C / PLM-mean-cell-C)
  set-current-plot-pen "PLM-Total"
  plotxy time ((PLM-cell-C + PLM-dead-cell-C) / PLM-mean-cell-C)
  
  ask turtles [set q-Sr q * Sr]
  let ave-q sum [q-Sr] of turtles with [alive = 1] / sum [Sr] of turtles with [alive = 1]
  ifelse count turtles with [alive = 0] > 0
  [set ave-q-dead sum [q-Sr] of turtles with [alive = 0] / sum [Sr] of turtles with [alive = 0]]
  [set ave-q-dead 0]
  
  set-current-plot "Quota"
  set-current-plot-pen "Alive"
  plotxy time ave-q * 10 ^ 3
  set-current-plot-pen "Dead"
  plotxy time ave-q-dead * 10 ^ 3
  set-current-plot-pen "PLM"
  plotxy time PLM-q * 10 ^ 3

  set sumS S * Vol
  ask turtles with [alive = 1] [set q-m-Sr q * m * Sr]
  ask turtles with [alive = 0] [set q-m-Sr q * m * Sr]
  set sumq (sum [q-m-Sr] of turtles with [alive = 1] + sum [q-m-Sr] of turtles with [alive = 0]) / Vol
  let S-t sumq + S  
  
  set-current-plot "Extracellular P"
  set-current-plot-pen "S"
  plotxy time S * 10 ^ 3
  set-current-plot-pen "S-PLM"
  plotxy time PLM-S * 10 ^ 3
end

to update-plot-CDF
set-current-plot "Internal P Content CDF"
clear-plot
set-current-plot-pen "Model"
let sumSr sum [Sr] of turtles
let cdfSr 0.0 
foreach sort-on [q * m] turtles
  [ ask ? [ 
      set cdfSr cdfSr + Sr
      set CDF-probability cdfSr / (sumSr + 1)
       
 plotxy (q * m * 10 ^ 15) (cdfSr / (sumSr + 1))
 ifelse alive = 1
[set-plot-pen-color 55]
[set-plot-pen-color 15]]]
 
 set-current-plot "Internal P Content CDF"
 set-current-plot-pen "Data"
 let number 0
 while [number < length Data-x] 
 [plotxy (item number Data-x) (item number Data-y) set number number + 1]
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Move Routines;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
to move-turtles-1
  ask turtles [
    right random 360
    forward random 10
    set color scale-color green (q * m) 0 8e-15
  ]
end

to move-turtles-2
  ask turtles [
    set xcor m * (max-pxcor  - min-pxcor) / (2 * m0 - m0) + 1
    set ycor q * (max-pycor  - min-pycor) / ((qmax-mean + (qmax-mean * qmax-CV * 2)) * 1.2  - q0)
     set color scale-color green q 0 .05
  ]
  ;ask patches [
  ;if-else LI = 1
  ;[
  ;  set pcolor 9.9
  ;]
  ;[
  ;  set pcolor 0
  ;]
  ;]
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Output;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
to output-data
  set file-number 0
  check-file-2
  file-open (word (filename-2))
  ask turtles [
     file-type q file-type "," file-type m file-type "," file-type Sr file-type "," file-type alive file-type ","
  ]
  file-close
end

to check-file
 set file-number file-number + 1
 set file-prefix "Chemostat_Model_A12_"
 set file-suffix ".csv"
 set filename (word file-prefix (file-number) file-suffix)
if (file-exists? (word (filename))) [check-file]
end

to check-file-2
 set file-number file-number + 1
 set file-prefix "Heterogeneity_Model_Full_Output_M2_"
 set file-suffix ".csv"
 set filename-2 (word file-prefix (file-number) file-suffix)
if (file-exists? (word (filename-2))) [check-file-2]
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Statistics;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
to calculated-stats
ask turtles [set q-w q * Sr]
set weighted-average-q sum [q-w] of turtles / sum [sr] of turtles
set weighted-average-m sum [m * Sr] of turtles / sum [sr] of turtles
set weighted-std-dev-q (sum [sr * (q - weighted-average-q) ^ 2] of turtles / (((count turtles - 1) * sum [sr] of turtles) / count turtles)) ^ .5
set weighted-std-dev-m (sum [sr * (m - weighted-average-m) ^ 2] of turtles / (((count turtles - 1) * sum [sr] of turtles) / count turtles)) ^ .5
set weighted-CV weighted-std-dev-q / weighted-average-q
set weighted-average-qm sum [q * m * Sr] of turtles / sum [sr] of turtles
set weighted-Var sum [Sr * (q * m - weighted-average-qm) ^ 2] of turtles / ((count turtles - 1) * sum [sr] of turtles) / count turtles
set weighted-Var-m sum [Sr * (m - weighted-average-m) ^ 2] of turtles / ((count turtles - 1) * sum [sr] of turtles) / count turtles
end

to Root-Mean-Squared-Error
  let number 0
  let temp 0
  while [number < length Data-x] 
  [;type item number Data-y type ","
    ask min-one-of turtles [abs(CDF-probability - item number Data-y)] 
    [set temp temp + (q * m * 10 ^ 15 - item number Data-x) ^ 2
      ;type CDF-probability type ","
      ;type item number Data-x type ","
      ;print q-m * 10 ^ 15
      ]
     set number number + 1
    ]
  set RMSE (temp / (number + 1)) ^ .5
end

to check-Variability
  ifelse Upmax-Variability
  [
    set uPmax-CV 0
    ]
  [
    set uPmax-CV 0
    ]
  ifelse q0-Variability
  [
    set q0-CV 0.15
    ]
  [
    set q0-CV 0
    ]
  ifelse uR-Variability
  [
    set uR-CV 0.2
    ]
  [
    set uR-CV 0
    ]
  ifelse m0-Variability
  [
    set m0-CV 0
    ]
  [
    set m0-CV 0
    ]
  ifelse uDbase-Variability
  [
    set uDbase-CV 0.1
    ]
  [
    set uDbase-CV 0
    ]
  ifelse qD-Variability
  [
    set qD-CV 0.15
    ]
  [
    set qD-CV 0
    ]
  ifelse VAmax-Variability
  [
    set VAmax-CV 0.1
    ]
  [
    set VAmax-CV 0
    ]
  ifelse Km-Variability
  [
    set Km-CV 0.05
    ]
  [
    set Km-CV 0
    ]
  ifelse qmax-Variability
  [
    set qmax-CV 0.6
    ]
  [
    set qmax-CV 0
    ]
  ifelse kw-Variability
  [
    set kw-CV 0
    ]
  [
    set kw-CV 0
    ]
  ifelse Sf-qm-Variability
  [
    set Sf-qm-CV 0.4
    ]
  [
    set Sf-qm-CV 0
    ]
  ifelse Sf-m-Variability
  [
    set Sf-m-CV 0.05
    ]
  [
    set Sf-m-CV 0
    ]
end

;;OPT::;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;OPT::;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;OPT::;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Optimization Routine;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;OPT::;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;OPT::;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;OPT::to optimize
;;OPT::  set filename "Heterogeneity_Model_A2_7.csv"
;;OPT::  set t-end-opt 15000
;;OPT::  set i 0
;;OPT::  set j 1
;;OPT::  set k 0
;;OPT::  set opt-step-size .5
;;OPT::  set opt-runs-step 3
;;OPT::  set shuffle-order shuffle shuffle-order
;;OPT::  set var-1 item 0 shuffle-order
;;OPT::  set var-2 item 1 shuffle-order
;;OPT::  set var-3 item 2 shuffle-order
;;OPT::  set var-4 item 3 shuffle-order
;;OPT::  set var-5 item 4 shuffle-order
;;OPT::  set var-6 item 5 shuffle-order
;;OPT::  set var-7 item 6 shuffle-order
;;OPT::  set var-8 item 7 shuffle-order
;;OPT::  set var-9 item 8 shuffle-order
;;OPT::  set var-10 item 9 shuffle-order
;;OPT::  set var-11 item 10 shuffle-order
;;OPT::  set var-12 item 11 shuffle-order
;;OPT::  let go-loop 0
;;OPT::  set opt-run-number 1
;;OPT::  while [k < opt-runs-step]
;;OPT::  [
;;OPT::  setup-opt
;;OPT::  set-variables
;;OPT::  set go-loop 0
;;OPT::  while [go-loop < t-end-opt] [opt-go set go-loop go-loop + 1]
;;OPT::  determine-rank
;;OPT::  Root-Mean-Squared-Error
;;OPT::  set k k + 1
;;OPT::  set opt-step-total opt-step-total + RMSE
;;OPT::  ]
;;OPT::  set RMSE opt-step-total / k
;;OPT::  set RMSE-min RMSE
;;OPT::  file-open (word (filename))
;;OPT::  file-type "0" file-type "," file-type RMSE file-type "," file-type RMSE-min file-type "," file-type item var-1 opt-fact file-type "," file-type item var-2 opt-fact file-type "," file-type item var-3 opt-fact file-type "," file-type item var-4 opt-fact file-type "," file-type item var-5 opt-fact file-type "," file-type item var-6 opt-fact file-type "," file-type item var-7 opt-fact file-type "," file-type item var-8 opt-fact file-type "," file-type item var-9 opt-fact file-type "," file-type item var-10 opt-fact file-type "," file-type item var-11 opt-fact file-type "," file-type item var-12 opt-fact file-type ","
;;OPT::  file-close
;;OPT::  while [j <= 80] [
;;OPT::    set var-current item i opt-fact
;;OPT::    ifelse random-float 1 <= .5
;;OPT::    [set opt-fact replace-item i opt-fact (var-current + opt-step-size)]
;;OPT::    [set opt-fact replace-item i opt-fact (var-current - opt-step-size)]
;;OPT::    if item i opt-fact < 0 [set opt-fact replace-item i opt-fact 0]
;;OPT::    set k 0
;;OPT::    set opt-step-total 0
;;OPT::    while [k < opt-runs-step]
;;OPT::    [
;;OPT::    setup-opt
;;OPT::    set-variables
;;OPT::    set go-loop 0
;;OPT::    while [go-loop <= t-end-opt] [opt-go set go-loop go-loop + 1]
;;OPT::    determine-rank
;;OPT::    Root-Mean-Squared-Error
;;OPT::    set k k + 1
;;OPT::    set opt-step-total opt-step-total + RMSE
;;OPT::    ]
;;OPT::    set RMSE opt-step-total / k
;;OPT::    ;update-plot-cdf
;;OPT::    file-open (word (filename))
;;OPT::    file-type opt-run-number file-type "," file-type RMSE file-type "," file-type RMSE-min file-type "," file-type item var-1 opt-fact file-type "," file-type item var-2 opt-fact file-type "," file-type item var-3 opt-fact file-type "," file-type item var-4 opt-fact file-type "," file-type item var-5 opt-fact file-type "," file-type item var-6 opt-fact file-type "," file-type item var-7 opt-fact file-type "," file-type item var-8 opt-fact file-type "," file-type item var-9 opt-fact file-type "," file-type item var-10 opt-fact file-type "," file-type item var-11 opt-fact file-type "," file-type item var-12 opt-fact file-type ","
;;OPT::    file-close
;;OPT::  ifelse RMSE < RMSE-min
;;OPT::    [
;;OPT::      set RMSE-min RMSE
;;OPT::    ]
;;OPT::      [
;;OPT::        set opt-fact replace-item i opt-fact var-current 
;;OPT::        ]
;;OPT::    set i i + 1
;;OPT::    
;;OPT::    set opt-run-number opt-run-number + 1
;;OPT::    if i > (length opt-fact - 1) 
;;OPT::    [
;;OPT::      set i 0 
;;OPT::      set j j + 1 
;;OPT::      ]
;;OPT::  ]   
;;OPT::end

;;OPT::to set-variables
;;OPT::  set uPmax-CV (uPmax-CV-base * item var-1 opt-fact)
;;OPT::  set q0-CV (q0-CV-base * item var-2 opt-fact)
;;OPT::  set uR-CV (uR-CV-base * item var-3 opt-fact)
;;OPT::  set m0-CV (m0-CV-base * item var-4 opt-fact)
;;OPT::  set uDbase-CV (uDbase-CV-base * item var-5 opt-fact) 
;;OPT::  set qD-CV (qD-CV-Base * item var-6 opt-fact)
;;OPT::  set VAmax-CV (VAmax-CV-base * item var-7 opt-fact)
;;OPT::  set Km-CV (Km-CV-base * item var-8 opt-fact)
;;OPT::  set qmax-CV (qmax-CV-base * item var-9 opt-fact)
;;OPT::  set Kw-CV (Kw-CV-base * item var-10 opt-fact)
;;OPT::  set Sf-qm-CV (Sf-qm-CV-base * item var-11 opt-fact)
;;OPT::  set Sf-m-CV (Sf-m-CV-base * item var-12 opt-fact)
;;OPT::end
;;OPT::
;;OPT::to setup-opt
;;OPT::  clear-turtles
;;OPT::  set LI 1
;;OPT::  set time 0
;;OPT::  set tlnext dtL
;;OPT::  set tdnext time + dt-d
;;OPT::  set S S0
;;OPT::  crt n-start
;;OPT::    [
;;OPT::    set shape "circle"
;;OPT::    set color green
;;OPT::    set size 1
;;OPT::    set xcor random-xcor
;;OPT::    set ycor random-ycor 
;;OPT::    set uPmax kreft uPmax-Mean uPmax-CV 0 99
;;OPT::    set q0 kreft q0-Mean q0-CV 0 99
;;OPT::    set uR kreft uR-Mean uR-CV 0 99
;;OPT::    set m0 kreft m0-Mean m0-CV 0 99
;;OPT::    set uDbase kreft uDbase-Mean uDbase-CV 0 99
;;OPT::    set qD kreft qD-Mean qD-CV 0 99
;;OPT::    set VAmax kreft VAmax-Mean VAmax-CV 0 99
;;OPT::    set Km kreft Km-Mean Km-CV 0 99
;;OPT::    set qmax kreft qmax-Mean qmax-CV 0 99
;;OPT::    set kw kreft kw-Mean Kw-CV 0 99
;;OPT::      
;;OPT::    set m m0-Mean
;;OPT::    set size m * m-plot-scale
;;OPT::    set q  q0-Mean * 5
;;OPT::    set Sr C0 * 1000 / n-start
;;OPT::    set age 0
;;OPT::    set alive 1 
;;OPT::    set h (((m / p) * 4) / (pi * d ^ 2))
;;OPT::    set A .5 * pi * d ^ 2 + pi * d * h
;;OPT::    set CDF-probability 0
;;OPT::    ]
;;OPT::end
;;OPT::
;;OPT::To opt-go
;;OPT::  tick
;;OPT::  ask turtles with [alive = 1] [uptake-rate]
;;OPT::  ask turtles [excretion-rate]
;;OPT::  ask turtles [mass-balance-cells]
;;OPT::  mass-balance-chemostat
;;OPT::  ask turtles with [alive = 1] [growth-rate]
;;OPT::  ask turtles with [alive = 1] [Grow]
;;OPT::  ask turtles with [alive = 1] [update-Age]
;;OPT::  ask turtles with [alive = 1] [Divide]
;;OPT::  ask turtles with [alive = 1] [mortality]
;;OPT::  if count turtles < 1 [stop]
;;OPT::  split-combine
;;OPT::  if time > tdnext
;;OPT::  [dilute]
;;OPT::  ;update-plot-CDF
;;OPT::  ;update-plot
;;OPT::  set time time + dt
;;OPT::  if time > tlnext [update-light]
;;OPT::end
;;OPT::
;;OPT::to determine-rank
;;OPT::let sumSr sum [Sr] of turtles
;;OPT::let cdfSr 0.0 
;;OPT::foreach sort-on [q * m] turtles
;;OPT::  [ ask ? [ 
;;OPT::      set cdfSr cdfSr + Sr
;;OPT::      set CDF-probability cdfSr / (sumSr + 1)
;;OPT::      ]
;;OPT::  ]
;;OPT::end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;PLM Routine;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
to PLM
   set PLM-V PLM-Vmax * PLM-S / (PLM-Km + PLM-S) * (PLM-qmax - PLM-q) / (PLM-qmax - PLM-q0)
   set PLM-W PLM-q * PLM-kw
   set PLM-uG PLM-uPmax * (1 - PLM-q0 / PLM-q) * LI - uR-Mean
   set PLM-uD PLM-uDbase * PLM-qd / PLM-q
   set PLM-q PLM-q + ((PLM-V - PLM-W) - PLM-uD * PLM-q - PLM-q * PLM-ug) * dt
   set PLM-dead-cell-C PLM-dead-cell-C + (PLM-uD * PLM-cell-C - Flow / Vol * PLM-dead-cell-C) * dt
   set PLM-cell-C PLM-cell-C + (PLM-cell-C * (PLM-uG - PLM-uD) - Flow / Vol * PLM-cell-C) * dt
   set PLM-S PLM-S + (Flow / Vol * S-in - Flow / Vol * S + (PLM-W - PLM-V) / .001 * PLM-cell-C) * dt
end