Grid Molecular Concentration
Often cells depend on the nutrient molecules that exist in the extracellular space to grow. In this model, we introduce nutrient uptake from the simulation domain. Modeling diffusion in extracellular space requires the use of built-in Partial Differential Equations solver in model_routine_grid.cpp. Additionally, we add cell death which depends on molecular concentration to model necrotic death.
This tutorial is organized in the following manner:
Section 2. Model Descriptions explains the 3 files changed from the previous tutorial:
Section 3. Model Output describes the expected output and animation of the model.
We add an additional enumeration, CELL_MODEL_REAL_ALIVE to cell_model_real_e to indicate whether a cell is dead or not. This will be used in updateSpAgentOutput() in model_routine_output.cpp to mark different colors for live and dead agents.
We want to model extracellular space where the entire simulation domain contains nutrients. It requires a number of new grid-specific model variables:
GRID_MOLECULAR_CONCENTRATION A_ALPHA A_NUTRIENT_VOLUME A_DOUBLING_TIME A_CELL_UPTAKE_RATE UPTAKE_PCT_INC_RATIO A_K diffusible_elem_e grid_model_real_e
GRID_MOLECULAR_CONCENTRATION is the concentration, phi, that will be initialized for every unit box in the simulation domain.
A_ALPHA is the amount of nutrient it takes for a cell to reach its maximum cell size and divide into mother and daughter cells. A_ALPHA is used with A_DOUBLING_TIME, the time it takes for a cell to double, to compute constant A_CELL_UPTAKE_RATE. The amount of absorbed nutrient mass is proportional to the uptake rate and uScale. Details on uptake rate and uScale will be further discussed in section, model_routine_grid.cpp.
A_NUTRIENT_VOLUME is the proportional mass of nutrient.
UPTAKE_PCT_INC_RATIO is used for gentle increase in uptake rate when cells are initially born. Initial “mature” cells have the CELL_MODEL_REAL_UPTAKE_PERCT of 1.0, and may has normal A_CELL_UPTAKE_RATE. However, when cells are immediately born, they should not have the same A_CELL_UPTAKE_RATE as fully mature cells, as they are not the same in size, and thus absorb nutrients at smaller rates. Hence, at cell division, mother and daughter cells are initialized CELL_MODEL_REAL_UPTAKE_PERCT of 5.0, and the uptake percentages are gradually increased by each time-step by the factor of UPTAKE_PCT_INC_RATIO. CELL_MODEL_REAL_UPTAKE_PERCT may not exceed 1.0.
A_K is the constant molecular concentration at which cells grow at half of their maximal growth rate.
In diffusible_elem_e, our diffusible nutrient of interest that exists in all simulation domain is enumerated as DIFFUSIBLE_ELEM0. model_routine_config.cpp
In solving for the diffusion of nutrients in the simulation domain, the Biocellion’s Partial Differential Equations setup, updatePDEInfo() has to be configured.
In updatePDEInfo(), create an ‘PDEInfo’ object and fill out associated members, pdeType, numLevels, numTimeSteps, callAdjustRHSTimeDependentLinear, mgSolveInfo, advectionInfo, and spllitingInfo. We do not use advection equation or splitting scheme in solving PDE’s, hence they are left as dumm variables.
In updatePDEInfo(), create a ‘GridPhiInfo’ object to set the element index (elemIdx), grid concentratation name (name), and boundary conditions at the six ends of the simulation domain (aa_bcType and aa_bcVal). We use the Neumann boundary condition at values 0, as we do not want any change in flux at any of the boundaries (we want to “insulate” the nutrients in the simulation domain).
In updateIfGridModelVarInfo(), define the three variable properties of each grid unit box: GRID_MODEL_REAL_COLONY_VOL_RATIO, GRID_MODEL_REAL_U_SCALE, and GRID_MODEL_REAL_RHS .
In updateFileOutputInfo(), set fileOutputInfo.particleOutput and fileOutputInfo.v_gridPhiOutput as ‘true’ in order for Biocellion to output PVTU and VTI files for visualization in Paraview.
We add a few grid model information in updateSummaryOutputInfo() to observe change in number of cells, and sum/average/min/max of the nutrient concentration during simulation output. However, it is up to the modelers to determine the model agent/grid model outputs. Less summary outputs at big or no intervals will result in faster model simulation.
From the previous addSpAgents() model script, add the additional cell state, CELL_MODEL_REAL_ALIVE with the initial value of 1.0 (alive).
Function updateSpAgentsState() is heavily altered from the previous tutorial models. The biomass of a cell is now dependent on the nutrient uptake of nutrients from the neighboring grid states. We consider the 9 unit boxes surrounding and including the the specific agent of interest. Loop over the 9 unit boxes (nested for loop) to find the “proportional” cumulative uScale. uScale is then used to find growthVolume. Growth of a cell at each BASELINE_TIME_STEP is proportional to the A_CELL_UPTAKE_RATE and uScale and uptakePct. Add computed growthVolume to oldVolume. Compute newRadius of the cell from newVolume, and also newBiomass using A_CELL_DENSITY.
uptakePct is used to tune for slower uptake (to reduce excessive diffusion) at cell birth. When cells are first born, they have the uptakePct of .5. uptakePct increases by the factor of UPTAKE_PCT_INC_RATIO at each BASELINE_TIME_STEP. However, uptakePct can not be greater than 1.0, so when uptakePct is less than 1.0, uptakePct = 1.0.
Update the the new agent states, newRadius and newBiomass.
In updateSpAgentBirthDeath(), set disappear = false, as we would like to visualize dead cells in Paraview.
In adjustSpAgent(), add additional way of cell death depending on remaining nutrient in the extracellular space. As cells consume more and more nutrients, and they become scarcer, cells will start dying. This models necrotic death that may resembles tumor growth:
(tmp2 < (1 - v_gridPhiNbrBox[DIFFUSIBLE_ELEM0].getVal(0,0,0)/GRID_MOLECULAR_CONCENTRATION)
Set agent states CELL_MODEL_REAL_UPTAKE_PERCT and CELL_MODEL_REAL_ALIVE as 0.0.
In divideSpAgent(), include the new cell state, CELL_MODEL_REAL_ALIVE with value of 1.0.
model_routine_grid.cpp handles all the functions that occur in extracellular space.
Similarly to addSpAgents(), initialize molecular concentrations and grid box properties in initIfGridVar(). This function is called for every single unit box in the simulation domain. The grid molecular concentration, v_gridPhi, is set with the value GRID_MOLECULAR_CONCENTRATION for every box. Set additional unit box grid properties, uScale and rhs to 0.0.
In updateIfGridVar(), we determine the diffusion and change in nutrient concentrations in the extracellular space that update at each NUM_STATE_AND_GRID_TIME_STEPS. The goal here is to compute the grid model consumption rate, rhs. Given the diffusion equation , the consumption rate is , and is defined . a_uPrime is cumulative A_CELL_UP_UPTAKE_RATE * uptakePct * ratio. a_uPrime is computed cumulatively (nested for loop) in the agent’s unit box and its 8 neighboring unit boxing, similarly to adjustSpAgent(). a_Uscale is the Monod’s equation: concentration of the limiting nutrient, phi, divided by phi + A_K, where A_K is the concentration at which half maximal growth is achieved. Save the updated states, GRID_MODEL_REAL_U_SCALE, GRID_MODEL_REAL_RHS, and GRID_MODEL_REAL_COLONY_VOL_RATIO.
Initialize updateIfGridRHSLinear() with GRID_MODEL_REAL_RHS. This function adds the consumption rate, , to the Diffusion equation.
In updateSpAgentOutput(), use CELL_MODEL_REAL_ALIVE to set different colors for dead and live cells in the .pvtu output.
Use updateSummaryVar(), to output user-set NUM_GRID_SUMMARY_REALS.
Load VTI and PVTU files in Paraview to simultaneously visualize growing/dying agents and molecular concentration, phi. Refer to the attached animation.