Full Text
Contents
1. CHECK QW GNLERWS Ze os 3 6 Setup fora Simulation Instance Simulating Microbial Community Patterning Using Biocellion 251 MODEL END return Code 5 Model routine to set the discrete agent color variable for visualization Biocellion asks users to provide specifics of a simulation instance e g simulation domain size output directory in an xml file See SBIOCELLION ROOT framework main yeast patterning 5um xml SBIOCELLION ROOT framework main yeast patterning 20um xml or SBIOCELLION ROOT framework main yeast pat terning 40um xml for examples for 5 20 or 40 um grid spacing respectively We set the required parameters first 1 Set the number of base line steps to execute We set this number to 60 000 500 h zzme step num baseline 60000 gt 2 Set the simulation domain size lt domain x 128 y 128 z 4864 in case we adopt 5 um interface grid spacing As we are using three AMR levels with the refinement ratio of 4 with 5 um interface grid spacing refinement ratio is also set in this xml file simulation domain size should be a multiple of 64 We set the domain size in the x and y directions slightly smaller than the size in 5 while setting the size in the z direc tion slightly larger than the size used in 5 Set the simulation initialization method and the partition size We set the initialization method to initialize within the code and set the partiti2. Simulating Microbial Community Patterning Using Biocellion 2493 higher performance See Note 1 Bzocellion requires the Intel Thread Building Blocks library freely available from the thread building blocks homepage http threadingbuildingblocks org l Download the most recent stable version of Intel Thread Building Blocks library version 4 2 or later is required to the target system 2 Unzip the downloaded tarball 3 Update the LD LIBRARY PATH Linux environment variable to include the TBB library directory in TBB 4 2 this is TBB ROOT lib intel64 gcc4 1 1 Unzip the Bzocellion tarball 2 Open Makefile common under the Bzocellion root directory 3 Update BIOCELLION ROOT to point to the Bzocellion root directory 4 Try make under the libmodel directory This should compile the model library SBIOCELLION ROOT libmodel interface libmodel so Biocellion users provide model specifics by filling in a set of C functions defined in five files under the SBIOCELLION ROOT libmodel model directory model routine config cpp model routine agent cpp model routine mech _ Intre6teopps model routine grid cpp and model routine output cpp These files include model routines to initialize the model update discrete agent states simulate direct physico mechanical interactions update the state of the extracel lular space and set simulation output respectively The entire model code for the Momeni et al s mo
3. VectorshRbEAh s amp V ogridPRLi const VecLtor REAL S v gridModelReal0 const Vector lt REAL gt amp v gridModelReall const Vector lt S32 gt b v grrdModellutO0 const Vector amp 552 2 Vv gridModelIntl1 REAL amp gridBeta MODEL START REAL z0 REAL vIdx0 2 0 5 gt IF GRID SPACING REAL zl REAL vIdx1 2 0 5 gt IF GRID SPACING REAL gridBetaO0 250 Seunghwa Kang et al REAL gridBetal TE zD S AGAR HEIGHT gridBeta0 A DIFFUSION COREF AGAR elemraoax else REAL scale REAL ubAgentData0 v spAgent size REAL UB FULL CEBL CNT if scale gt 1 0 scale 1707 gridBetaU A DIFFUSION COEFF COLONY Lelemidx scale dE d ue AGAR HEIGHT gridBetal A DIFFUSION CORFF AGAR elemrdax else REAL scale REAL ubAgentDatal v spAgent size REAL UB PUBL CELL CNTs if scale gt 1 0 scale 1 0 GridBeta r DLERUSION COREF COLONY Lelemidx scale if gridBetaO gt 0 0 amp amp gridBetal gt 0 0 gridBeta 1 0 1 0 gridBeta0 1 0 grid Betaly wco0u5 y harmonic mean else gridBeta 0 0 MODEL END return Code 4 Model routine to set diffusion coefficient between two adjacent grid boxes in the interface grid void ModelRoutine updateSpAgentOutput const VIdx amp vldx const SpAgent amp SpAgent REAL amp Color Vector lt REAL gt amp V Xetra 4 MODEL START color spAgent state getType
4. and estimated diffusion rates does not exceed the lysine con centration of the grid box Model routines setting PDE parameters and model routines updating individual cell states can access these values to set nutrient uptake rates this is necessary to assure that the total amount of nutrients consumed in solving PDEs coincide with the amount consumed by cells in updating individual cell states We save 7rbs and rhs 4 to avoid computing the rates again when setting PDE parameters A G strain cell releases lysine on death A relatively large amount of lysine is released in a short amount of time and this forms a steep concentration gradient followed by a rapid gradient change due to diffusion Accurately computing this transient gradient change requires a small time step size However cells react to concentration changes only gradually so accurately computing the transient gradi ent change has little impact on simulation output We have decided to spread released lysine to six neighboring boxes in the x y and z directions to lower the initial concentration gradient We implement this by updating the 775 variable of a neighboring grid boxes Figure 3 shows a simulation output Simulation time is highly dependent on grid resolution and with 40 um interface grid spac ing comparable to 50 mum grid spacing used in the Momeni et al s work 5 to simulate 500 h of cell growth a single simula tion run takes 6 5 h on a workstation
5. simulation domain in the z direction updatePDEBufferBdryType Set the boundary type between an interface grid partition and a PDE buffer grid partition The current version of Bzocellion provides only one option PDE BUFFER BDRY TYPE HARD WALL and discrete agents are not allowed to pass the boundary This function is irrelevant to PDE boundary conditions updateTimeStepInfo Update time step sizes Set the baseline time step size to 30 s and split a single baseline time step to 30 state and grid time steps Interested readers can experiment with different time step sizes updateSyncMethod Update synchronization methods when a single variable is updated by multiple model routine calls see Note 2 The only kind of short range cell cell direct mechanical interaction we consider is cell cell shoving so the synchronization method for extra mechanical interactions is irrelevant Set the synchronization method for grid variable updates to SYNC METHOD DELTA see Note 2 updateSpAgentInfo Set discrete agent types We consider three cell types R and G strain cells and dead cells dead cells do not grow and divide A dead cell has one model specific variable storing the amount of lysine in the cell This variable 1s used to calculate the amount of lysine to be released into the extracel lular space updatePDEInfo Set the grid state variables for lysine and adenine concentrations that are updated by solving PDEs An AM
6. zpdatelfGrid R HSLinear sets the reaction term based on the lysine and adenine production and consumption rates for a grid box in the interface grid No cells reside in the PDE buffer region and updatePDEBufferR HSLin ear sets the reaction term to 0 0 7 updatelfGridAMRTags In solving PDEs Bzocellion users can apply different grid spacings for different regions the finest grid spacing is the interface grid spacing Users tag each box in the interface region with a desired AMR level Boxes in the PDE buffer region are assumed to be tagged with the coarsest level Bzocellion using CHOMBO 7 generates an AMR hierarchy based on this information We tag the boxes contain ing cells and the top 40 um in the z direction in the agarose cylinder with the finest level We tag the remaining boxes with the coarsest level 3 5 Simulation Model routines controlling simulation outputs are defined in Output model routine output cpp l updateSpAgentOutput Color discrete agents We color each discrete agent based on the cell type We do not need to update extra output variables as we are mapping a discrete agent to a sphere See Code 5 for our implementation for the Bzocellion framework void ModelRoutine updateIfGridBetaInIfRegion Const 532 elLemldx Const 932 dim onset ViCcxe vLaxU const Vidxe YVIdOdXl const UBAgentData amp ubAgentData0 const UBAgentData amp ubAgent Datal const Vector lt REAL gt amp V COgIIGOPRhEU Conse
7. in slow turn around time In contrast with efficient parallel software a single simulation run can be partitioned across multiple computing cores dramatically reducing the turn around time Biocellion is a high performance software framework that enables computational biologists without parallel computing expertise to exploit the power of parallel computers with only moderate programming effort We briefly describe Momeni et al s model 5 and Bzocellion 6 and illustrate the necessary steps to adapt Momeni et al s model to Bzocellion as a case study Momeni et al 5 implemented a mathematical model with a range of parameters corresponding to different yeast strains engag ing in various interactions We consider only one instance in this article Modification for other strains and interaction types is straightforward Two strains of yeast grow on top of a 24 mm thick agarose cylinder One strain consumes lysine and also secretes adenine at a constant rate say the R strain The other strain consumes adenine and releases lysine on death say the G strain The two strains grow cooperatively and Momeni et al demonstrated that strong cooperation promotes intermixing of the two strains in a Simulating Microbial Community Patterning Using Biocellion 235 three dimensional community We list model specifics necessary to adapt the model to Bzocellton 1 Yeast cells grow on top of a 24 000 pm thick agarose cylinder Partial different
8. initPDEBufferPin Initialize grid state vari ables for the interface grid zmztIfGrid Var and the PDE buffer grid zzitPDEBufferPhi initIfGrid Var is called once per grid box in the interface grid and inttPDEBufferPhi is called once per grid box in the PDE buffer grid 2 updatelfGrid Var Update interface grid state variables based on model specific rules see Subheading 1 3 This function is called once per grid box in the interface grid void ModelRoutine computeForceSpAgent const VIdx amp vidx0 const SpAgent amp SpAgentO0 const VIdx amp vidxl const SpAgent amp spAgentl const VReal amp dir unit direction vector from spAgentl to spA gentO const REALA dist VReal amp force force on spAgent0 due to inter action with spAgentl force 248 Seunghwa Kang et al on spAgentl due to interaction with spAgentO0 has the same magnitude but the opposite direction if force has the same direction with dir two cells push each other if has the opposite direction two cells pull each Orne AAR 4 MODEL START REAL R spAgent0 state getRadius spA gentl state getRadius REA magy cor repulsive stores iof adhesive force it dist lt R u shoving to remove the overlap mag eee ec arse y j else adhesion mag 0 0 no adhesion j FOr S32 dim OF dim lt DIMENSION dimer 4 force dim mag dir dim j MODEL END TOES Code 3 Model routine t
9. Chapter 16 Simulating Microbial Community Patterning Using Biocellion Seunghwa Kang Simon Kahan and Babak Momeni Abstract Mathematical modeling and computer simulation are important tools for understanding complex interactions between cells and their biotic and abiotic environment similarities and differences between modeled and observed behavior provide the basis for hypothesis formation Momeni et al Elife 2 e00230 2013 investigated pattern formation in communities of yeast strains engaging in different types of ecological interactions comparing the predictions of mathematical modeling and simulation to actual patterns observed in wet lab experiments However simulations of millions of cells in a three dimensional com munity are extremely time consuming One simulation run in MATLAB may take a week or longer inhib iting exploration of the vast space of parameter combinations and assumptions Improving the speed scale and accuracy of such simulations facilitates hypothesis formation and expedites discovery Bzocellion is a high performance software framework for accelerating discrete agent based simulation of biological systems with millions to trillions of cells Simulations of comparable scale and accuracy to those taking a week of computer time using MATLAB require just hours using Bzocellion on a multicore workstation Bzocellion further accelerates large scale high resolution simulations using cluster computers by partitioning t
10. I The model library links to the Bzocellion core framework at runtime Biocellion output files can be visualized using Paraview http Www paraview org Biocellion has three computational modules to simulate 1 individual discrete agent behavior 2 direct physico mechanical interactions between discrete agents and 3 changes in the extra cellular environment Bzocellton imposes a grid on the simulation domain to represent the state of the extracellular environment In Momeni et al s model 5 cells reside on top of the 24 mm Simulating Microbial Community Patterning Using Biocellion 23 thick agarose cylinder and the region where yeast cells grow is a small fraction of the entire simulation domain The agarose cylinder is relevant only in tracking molecular concentrations in the model Maintaining data structures for all three computational modules for the entire agarose cylinder can waste a significant amount of computing and memory Bzocellion imposes two different types of grids to different parts of the simulation domain to avoid such waste Bzocellion imposes an interface grid on a region where all three computational modules are executed and computational modules communicate through this interface grid Bzocellion imposes a coarser PDE buffer grid on the region relevant only in solving PDEs e g tracking nutrient concentrations in the aga rose cylinder Biocellion decomposes the simulation domain into multiple par ti
11. Math Biol 58 4 5 657 687 2 Colella P Graves DT Johnson JN Johansen HS Keen ND Ligocki TJ Martin DF McCorquodale PW Modiano D Schwartz PO Sternberg TD Van Straalen B 2012 Chombo software package for AMR applications design document Lawrence Berkeley National Laboratory Berkeley CA 3 Ferrer J Prats C L pez D 2008 Individual based modelling an essential tool for microbiol ogy J Biol Phys 34 1 2 19 37 4 Galle J Loeffler M Drasdo D 2005 Modeling the effect of deregulated proliferation and apoptosis on the growth dynamics of epithelial cell populations in vitro Biophys J 88 62 75 Momeni B Brileya KA Fields MW Shou W 2013 Strong inter population cooperation leads to partner intermixing in microbial com munities Elife 2 e00230 Pacific Northwest National Laboratory 2013 Biocellion 1 0 User Manual 1 0 edition Accessed Jul 2013 Xavier JB Picioreanu C van Loosdrecht MCM 2005 A framework for multidimensional modelling of activity and structure of multispe cies biofilms Environ Microbiol 7 8 1085 1103
12. R scheme is used with three levels 5 um grid spacing or two levels 20 um grid spacing or 40 pm grid spacing The finest level has the grid spacing equal to the interface grid The size of a single PDE time step is set identical to the size of a single state and grid time step A zero flux boundary condition is applied in the z direction Boundary conditions in the x and y directions are irrelevant because periodic boundary conditions are imposed in zpdateDomainBdryType updatelfGrid ModelVarInfo Set extra model specific variables associated with the interface grid We add four variables per grid box see Subheading 1 3 updateR NGInfo Set one random number generator to get random numbers with the uniform distribution setPDEBuffer Set a partition as a PDE buffer partition if the top surface in the z direction of the partition is lower than the agarose cylinder height minus a small buffer region to 246 3 2 Individual Agent Behavior l 2 3 3 3 Physico Mechanical Interaction Between Agents Seunghwa Kang et al 12 exclude the top part of the agarose cylinder where molecular concentrations are highly localized setPDEBuffer is called once per partition setHalitable Set a grid box in the interface grid as a habitable box or an uninhabitable box We set the boxes in the agarose cylinder to be uninhabitable and cells are not allowed to move into an uninhabitable box Model routines related to simulati
13. bial Community Patterning Using Biocellion 239 production and consumption rates of extracellular molecules and this drives molecular concentration changes in the extracellular space These two modules can be coupled more tightly by splitting a single baseline time step into multiple state and grid time steps Variables associated with the grid imposed on the extracellular space can be updated either by model specific rules or by solving PDEs Users can update the variables by model specific rules at the beginning and at the end of each state and grid time step A single state and grid time step can be further partitioned to smaller PDE time steps to advance PDEs The original yeast patterning model 5 partitions the simulation domain on top of the agarose cylinder into a set of fixed size cubic boxes a box width is 5 um and a cell takes a single box If a cell divides the new cell tries to occupy one of the nearest neighboring boxes in the same z plane if there is an empty box within the confinement neighborhood of 5 cell radius The existence of con finement neighborhood was observed experimentally If there 1s no empty box within the confinement neighborhood in the same z plane the new cell occupies the box right on top of the mother cell box and all the other cell boxes on top of the mother cell box are pushed upward Using Bzocellion we represent each cell by a sphere Bzocellton allows users to map a discrete agent to a different s
14. cylinder surface at the beginning of the simulation Initial cell volume is set to a random value between one half of the maximum cell volume and the maximum cell volume to represent the range between a new born daughter and a fully grown cell The initial cell density is 500 cells per mm 8 Initial lysine and adenine concentrations are set to zero 9 Diffusion coefficients are 300 pm s in the agarose cylinder and 20 um s inside a yeast colony for both lysine and adenine according to experimental measurements Diffusion coeffi cients for grid boxes containing yeast cells are scaled down based on the total volume of the cells in a grid box a grid box with low cell volume a box in the colony air boundary has smaller diffusion coefficients than a box with high cell volume as implemented in the original model Biocellion s design goal is to accelerate a wide range of discrete agent based mathematical models of multicellular biological sys tems This is challenging because mathematical models of biologi cal systems vary significantly Bzocellton s approach is to separate model specifics from common computational and parallel pro gramming challenges Bzocellton asks users to provide model specif ics and Bzocellion handles the remaining computational and programming challenges Model specifics are expressed by the developer through modification of a library of C functions that comprise Bzocellion s Application Programming Interface AP
15. del 5 is available under SBIOCELLION ROOT libmodel model yeast patterning for interested readers Below we present examples of how a model is specified Model routines related to model configuration are defined in model POULIDne COntig cpp l updatelIfGridSpacing Set the interface grid spacing by filling the updatelfGridSpacing function body surrounded by MODEL START and MODEL END Code 1 The interface grid spacing should be equal to or larger than the maximum direct physico mechanical interaction distance We map a cell to a sphere shaped discrete agent and only consider cell shoving in evaluating short range mechanical interactions 244 Seunghwa Kang et al 20 um Agarose Gel Agarose Gel Fig 4 Comparing 5 um grid spacing and 20 um grid spacing Grid boxes outside the agarose cylinder with no cells have 0 diffusion coefficient If cell A secretes a metabolite consumed by cell B the secreted metabolite is delivered to cell B via diffusion through the agarose cylinder However if cell A and cell B are located in the same box with 20 um grid spacing cell B can directly consume the molecules secreted by cell A Two overlaping spheres are pushed apart to remove the overlap and the maximum mechanical interaction distance cannot exceed the maximum cell diameter 5 pum We start with the smallest interface grid spacing 5 um to minimize simulation artifacts see Fig 4 for an example Larger values reduce exec
16. e say Alysine is V x l4 Alysine where V is the Oy minimum cell volume the volume of a cell right after cell division or one half of the maximum cell volume o the amount of lysine required to produce a daughter cell is 2 fmol Similarly the volume of a G strain cell after consuming a certain Aadeni amount of adenine say Aadenine is V x o adenine Aa QC the amount of adenine required to produce a daughter cell is 1 fmol When a cell grows above the maximum cell volume the cell divides into two equal volume cells 5 For an R strain cell lysine uptake rate is v TL x T K A QA Pa For a G strain cell adenine uptake rate is v x T K s 236 Seunghwa Kang et al 1 2 Biocellion Overview The minimal population doubling times for lysine requiring and adenine requiring cells T and T4 respectively are 1 76 and 1 98 h respectively The Monod s constant Ki i e the concentration of lysine at which lysine requiring cells grow at their half maximal growth rate is 1 uM The Monad s constant KA for adenine requiring cells is 0 1 uM 4 and 4 are lysine and adenine concentrations in the extracellular space respec tively An R strain cell secretes 0 08 fmol of adenine per hour and a G strain cell releases 12 fmol of lysine on death 6 The death rates are 0 054 h and 0 018 h for an R strain cell and a G strain cell respectively 7 Cells are randomly distributed on the agarose
17. ete agents see Code 3 This model routine 1s called once per cell pair that is within the maximum direct mechanical interac tion distance at every baseline time step Force on two interacting cells 1s set based on the overlap between the two cells to remove the overlap by pushing the two cells apart void ModelRoutine updateSpAgentBirthDeath c onst VIdx amp vldx const SpAgent amp spAgent const AgentMechIntrctData amp mechIn trerData const Vector NbrtrBox 3 4 State Changes in the Extracellular Space Simulating Microbial Community Patterning Using Biocellion 247 REAL gt amp v gridPhiNbrBox elemIdx const Vector lt NbrBox lt REAL gt amp v gridModelRealNbrBox elemIdx const Vector NbrBox S32 gt amp v gridModelIntNbrBox elemIdx BOOL amp divide BOOL amp disappear MODEL START divide false disappear false if spAgent state getlType AGENT TYPE R CELL spAgent state getType AGENT TYPE G CELL R or G straincells if spAgent state getRadius gt MAX CELL RADIUS divide true j else dead cell CHECK SpAgent state getrype AGENT _ TYPE D CELL MODEL END return Code 2 Model routine to set whether a cell will divide disappear or neither divide nor disappear Model routines related to simulating state changes in the extracel lular space are defined in model routine grid cpp l znitIf Grid Var and
18. hape A sphere can be located anywhere in the simulation domain above the agarose cylinder and its radius changes to model cell growth When a cell grows just enough to overlap with another the model immediately introduces a force to push all spheres apart thus modeling cell shoving in packed regions Figure 2 shows that the concentration of adenine changes smoothly in the simulation domain except for the cell growth region and the top 40 um part of the agarose cylinder adenine concentration changes smoothly even just 40 um below the agarose cylinder top surface We set the partitions at the agarose cylinder region except for the top 40 um part as PDE buffer Bzocellion supports AMR which applies different grid resolutions to different parts of the simulation domain In generating an AMR hierarchy the region occupied by yeast cells and the top 40 um of the agarose cylinder are tagged with the finest grid spacing which is equal to the interface grid spacing A coarser grid is imposed on the air region and the bottom part of the agarose cylinder We set the baseline time step size to 30 s and split a single base line time step to 30 state and grid time steps to tightly couple cell metabolic rate change and nutrient concentration change in the extracellular space PDE time step sizes to advance PDEs updating lysine and adenine concentrations are set identical to the state and grid time step size 240 Seunghwa Kang et al The boundary be
19. he work to run on multiple compute nodes Bzocellion targets computational biologists who have mathemati cal modeling backgrounds and basic C programming skills This chapter describes the necessary steps to adapt the original Momeni et al s model to the Bzocellion framework as a case study Key words Discrete agent based modeling Partial differential equation Adaptive mesh refinement High performance computing Cell system simulation 1 Introduction Discrete agent based modeling maps a multicellular biological system to a collection of discrete agents Discrete agent based modeling has been widely used to model various biological systems 1 3 5 7 Momeni et al 5 studied spatial patterning in a com munity of yeast strains engaging in metabolic interactions through a combination of discrete agent based mathematical modeling and wet lab experiments Finding solutions using mathematical models often requires implementation of these models as a computer program Lianhong Sun and Wenying Shou eds Engineering and Analyzing Multicellular Systems Methods and Protocols Methods in Molecular Biology vol 1151 DOI 10 1007 978 1 4939 0554 6_16 Springer Science Business Media New York 2014 233 234 Seunghwa Kang et al 1 1 Yeast Patterning Model Description Producing a high performance implementation that anticipates and accommodates easy revisions and refinements of a model as it evolves is a time consuming task Computati
20. ial equations PDE are commonly used to model spatio temporal variation of molecular concentrations in the extracellular space We adopt this approach and PDEs model adenine and lysine concentration changes in the agarose cylinder and the yeast cell community Because the initial spa tial distribution of the two yeast strains 1s uniform a small frac tion of the plate area in wet lab experiment is already representative of patterns observed for the entire plate Thus only a fraction of the plate area in wet lab experiments is con sidered in simulation In Momeni et al s work 5 simulation domain height is set to include the entire agarose cylinder in the simulation area plus up to 300 um above the agarose cylin der where yeast cells grow 2 Periodic boundary conditions both for PDEs and cell move ments are assumed in the x and y directions A zero flux boundary condition is applied at the bottom end of the aga rose cylinder and at the top of the simulation domain mole cules and cells cannot cross the top and bottom planes of the simulation domain 3 We map a cell to a sphere instead of a fixed size cubic box in the original model using Bzocellzon The maximum cell diameter 3 is 5 um The maximum cell volume is LLL um Yeast cells push against other cells when they are packed together 4 We increase the volume of a sphere to model cell growth The volume of an R strain cell after consuming a certain amount of lysin
21. ion asks users to set the synchronization method to properly update grid variables when a single variable is updated by multiple model routine calls We set the synchronization method to SYNC_METHOD DELTA to set the value by summing the differences from the initial value when updated by multiple model routine calls e g if three different model routines set the value of variable rhsj to 3 5 and 9 respectively then 7 5 1s set to 3 549 assuming that the initial value is 0 We need to reset the sum of lysine uptake and secretion rates 755 to 0 0 to ensure that this scheme properly works We configure Bzocellion to invoke model routines to edit grid variables at the end of a state and grid time step to reset the variable to 0 0 Simulating Microbial Community Patterning Using Biocellion 253 Acknowledgements Support for this research was provided by the Extreme Scale Computing Initiative and the Fundamental and Computational Sciences Directorate as part of the Laboratory Directed Research and Development Program at Pacific Northwest National Laboratory PNNL Portions of this work were conducted using PNNL Institutional Computing at PNNL PNNL is operated by Battelle for DOE under contract DE ACO5 76RLO 1830 B M is a Gordon and Betty Moore Foundation fellow of the Life Sciences Research Foundation References 1 Byrne H Drasdo D 2009 Individual based and continuum models of growing cell populations a comparison J
22. n have more fine boxes than the user input to improve efficiency processing a large number of small boxes is inefficient and guarantee correctness coarsening a fine grid first and refining the coarsened fine grid should produce the original fine grid or see the proper nesting condition in ref 2 Figure 1 depicts the simulation domain left and the gener ated AMR hierarchy right in our experiment assuming 40 um grid spacing two AMR levels and the refinement ratio of 4 we simulate with 5 20 and 40 um grid spacings Multicellular biological system simulation combines multiple biological processes such as cell movement diffusion of molecules and cell metabolic rate change Different biological processes have 238 Seunghwa Kang et al The cell growth region 320 um The cell growth Interface 7 region with yeast cells 4 D W id Partition Sa EIZIZEZEZIZIZS Y YY gri EN ai p A40umgrid SEPT R spacing The top 40 um of the 7 YH The agarose cylinder 24 mm PDE buffer grid 160 um grid spacing Fig 1 Biocellion imposes two different types of grids to the simulation domain eff An interface grid is imposed on the partition covering the cell growth region to simulate both cells and the environment The remaining partitions in the agarose cylinder region are set as PDE buffer these partitions are relevant only in solving PDEs to track lysine and adenine concentrations U
23. neration ar refine ratio 4fill_ratio 0 5 gt See the Bzocellion user manual 6 for additional details Set the optional parameters affecting the accuracy of the mul tigrid method used in solving PDEs sj parabolic solve mg num pre 3 mg num post 3 mg num bottom 3 ji v or w y mg max ite rations 50 mg epsilon 12 mg_hang 8 mg norm threshold 20 gt See the Biocellion user manual 6 for additional details Biocellion users can set the framework to perform checks on model routine outputs or input arguments of Bzocellton utility functions called inside model routines by enabling CHECK FLAG DENABLE CHECK 1in BIOCELLION ROOT Makefile model This often allows users to easily identify model program bugs such as using a random number genera tor without initialization or accessing a C STL vector array variable outside the array length We recommend Bzocellion users enable this option to verify their model routines prior to running full simulations Once the model routines are verified users should disable this check to expedite simulation Lysine released by a dying G strain cell 1s spread to seven grid boxes The sum of lysine uptake and secretion rates for a grid box can be updated by seven different model routine calls a model routine to update grid state variables is invoked once for every grid box in the interface grid in a single round Biocell
24. ng individual agent behavior are defined in model routine agent cpp l addSpAgents Randomly spread cells on top of the agarose cylinder to initialize the simulation addSpAgents is called once per interface grid partition updateSpAgentState Increase cell size based on the nutrient consumption rates We use the scaled 9 values o k Subheading 1 3 to compute the nutrient consumption rates This function is called once for every discrete agent in every state and grid time step updateSpAgentBirthDeath divideSpAgent and adjustSpAgent Set whether a discrete agent will divide or disappear updateSpA JentBirthDeath see Code 2 updateSpAgentBirthDeath is called once for every discrete agent in every baseline time step A cell divides if its size exceeds the maximum cell size If a cell is set to divide divideSpAgent is called A cell divides into two cells in a random direction and the two resulting cells volume is one half of the original cell s volume in the ported model Dead cells remain in the simulation domain so no discrete agent is set to disappear If neither divide nor disappear is set adjustSpAgent is called and this model routine updates the cell displacement based on the sum of the forces on the cell Model routines related to simulating physico mechanical interac tions between discrete agents are defined in model routine mech IH EO WepD computeForceSpAgent Compute forces between pairs of dis cr
25. o compute force between two interacting discrete agents 3 updatelfGridKappa and updatePDEBufferKappa Set PDE parameter x for the interface grid and the PDE buffer grid for each grid box x represents the cell volume exclusion in diffu sion As we are not considering cell volume exclusion x is set to 1 0 0 volume exclusion 4 updatelfGrid Alpha and updatePDEBufferAlpha Set PDE parameter a a sets the decay rate We ignore lysine and ade nine decay and set a to 0 0 5 updatelfGridBetaInIfR egion updateIfGrid Beta PDE BufferBdry updatelfGridBetaDomainBdry updatePDE BufferBetaInPDEBuffer Region and updatePDEBufferBetaDomainBdry Set PDE parameter f p sets the diffusion coefficient f is set between two grid boxes sharing a face updateIfGrid BetaInIf Region and update PDE BufferBetaInPDEBuffer Region Different model routines are called at the boundary between the interface grid and the PDE buffer grid updateIfGrid Beta PDEBufferBdry and the simula tion domain boundary updatelfGridBetaDomainBdry and updatePDEBufferBetaDomainBdry Code 4 sets the diffusion coefficient between two adjacent grid boxes in the interface Simulating Microbial Community Patterning Using Biocellion 249 grid based on the yeast patterning model specifics by filling the function body of the predefined Bzocellion model routine updateIfGrid BetaInIfR egion 6 updatelfGridR HSLinear and updatePDEBuffer R HSLinear Set the PDE reaction term
26. on size to 64 init data partition size 64 src code gt Alternatively users can start from check point data 4 Set the output directory path interval file formats and the number of extra output variables for each discrete agent lt out put path output directory path interval 120 particle pvtu num extra 0 grid vim gt We also set several optional parameters relevant to executing the yeast patterning model on multicore PCs 1 Set standard output verbosity from 0 to 5 to 1 lt stdout ver bosity 17 2 Set the number of threads to 12 for a multicore PC with 12 hardware threads lt system num_node_groups 1 num_nodes_ per group 1 num sockets per node 1 max_load_imbal ance 1 2 num thre ads 12 gt num node groups num nodes per group num sockets per node and max load imbalance are irrelevant for multicore PCs and are ignored 252 Seunghwa Kang et al 4 Notes 3 Set the summary report Bzocellion allows users to print the summary of the interface grid variables see the Bzocellion user manual 6 for additional details AMR regridding and check point intervals lt interval summary 10 load balance 120 regridding 120 checkpoint 6007 load_balance is irrelevant to multicore PCs and is ignored Set the refinement ratio in applying AMR and other parame ters controlling AMR hierarchy ge
27. onal biologists often favor low programming effort over efficiency flexibility and even accuracy Therefore simulations tend to run much longer than performance optimized code Even incremental model updates can require significant code revision For example Momeni et al s model 5 partitions the simulation domain into fixed size cubic boxes and maps each cell to a single box Once a cell divides the daughter cell occupies one of the nearest neighboring boxes instantly pushing surrounding cells outward to free up the space This approximation of cell growth and shoving may be sufficient for the studied problem but it may not work well for more com plex problems e g this approach cannot model a cell division that produces two cells differing in size Parallel computers ranging from multicore PCs to leadership class supercomputers provide significantly larger computing capac ity than a single compute core This computing capacity can address the computational challenges in simulating complex models when harnessed by simulation software that serves widely varying multi cellular biological system modeling requirements and runs efficiently on parallel computers However implementing such software 1s dif ficult and time consuming Even if computational biologists have access to cluster computers without parallel computing the advan tage is limited to running multiple simulations in parallel one simulation per compute core resulting
28. ration in the extracellular space and K the con centration of metabolite at which half maximal consumption rate is achieved is 1 0 and 0 1 uM for lysine and adeneine respec tively Uscat sino and Us adenine Store the 9 values for lysine and o k adenine respectively We want to limit the total amount of lysine or adenine consumed by cells in a grid box to be lower than the amount of lysine or adenine in the box plus an estimation of the amount of lysine or adenine diffuse into the box within a single state and grid time step This prevents from becoming negative 2 Materials Simulating Microbial Community Patterning Using Biocellion 241 especially when is small without using a tiny time step size This approach is accurate as long as our estimation of the diffusion rate is accurate We use the explicit Euler method to estimate the amount of diffusion which gives a reasonably accurate estimation for our choice of the state and grid time step size 1 s the molec ular Concentration gradient does not change significantly within l s We scale to limit nutrient consumption If the sum of o K Spine Or Ths and the estimated diffusion rate multiplied by the state and grid time step size exceeds the lysine or adenine concentration of the box we reduce the consumption rate and scale down U sine OL Uscateadenine SO the net decrease of the lysine or adenine concentration based on the production consump tion
29. sers tag interface gird boxes with the desired AMR level We tag the interface grid boxes containing yeast cells the ight blue boxes and the interface grid boxes at the top 40 um of the agarose cylinder the ight red boxes with the finer AMR level The remaining interface grid boxes are tagged with the coarser level Biocellion generates an AMR hierarchy right based on this information Note that the generated AMR hierarchy has more fine boxes than the user input to satisfy the proper nesting condition 7 different time step requirements to simulate the processes with sufficient accuracy For example simulating diffusion of molecules by solving PDEs often requires a significantly smaller time step size than the time step size required to simulate cell movement To accommodate multiple time step size requirements in multicel lular biological system simulation Bzocellion uses multiple time step sizes to simulate different computational modules and to communicate across the modules The baseline time step is the larg est time step used to simulate direct physico mechanical interac tions and discrete agent birth death and movement The module computing direct physico mechanical interactions communicates with the other two modules once per baseline time step Discrete agent states and the state of the extracellular environment affect each other For example cell metabolic rate change affects the 1 3 Porting Overview Simulating Micro
30. tions users set the partition size Users can impose either an inter face grid or a PDE buffer grid for each partition When users wish to run Bzocellion on a cluster computer with multiple compute nodes each node has multiple compute cores Bzocellton creates multiple compute processes and every compute process works on a different set of partitions note that a single process can exploit multiple compute cores in a single compute node to work on a single parti tion Separate output files are created for different partitions Biocellion supports adaptive mesh refinement AMR to solve PDEs AMR generates multiple levels of grids with different grid spacings based on the spatial resolution requirements of different simulation domain subregions Bzocellton asks users to set the num ber of AMR levels and the refinement ratio between two consecu tive AMR levels if the refinement ratio is set to 4 the coarser level grid spacing is four times larger than the finer level grid spacing The finest grid spacing coincides with the interface grid spacing Biocellion users set the interface grid spacing The PDE buffer grid spacing equals the coarsest grid spacing in the AMR hierarchy Users tag interface grid boxes with the desired AMR level PDE buffer grid boxes are automatically tagged with the coarsest AMR level Bzocellion generates an AMR hierarchy which is used to solve PDEs based on this information Note that the generated AMR hierarchy ca
31. tween the cell growth region and the agarose cylinder P 3e 17 2 EJ 2 5e 17 z S 0 2000 4000 6000 8000 10000 12000 14000 16000 18000 20000 22000 24000 z unit um 1 73242e 16 5 4B424e 17 4 69383e 17 T 6e 16 d 666 17 5256 17 1 278 16 4542 17 56 17 Be 17 PELA 46 17 4566 17 ASe 17 45 a 65e 17 4 38168e 17 lle 17 L3 Fig 2 Adenine concentrations in the simulation domain unit umol um 1 uM 107 umol um The top figure depicts the adenine concentration along the z axis passing the center of the simulation domain with 0 and 24 000 being the bottom and the top of the agarose cylinder respectively The bottom figures show the ade nine concentration at the z normal planes right on top of the agarose cylinder right right below the agarose cylinder top surface center and 40 uum below the agarose cylinder top surface right respectively The adenine concentration has higher spatial variation near the agarose cylinder top surface The spatial variation of the adenine concentration is significantly lower even just 40 uum below the top surface We define four model specific variables say vhsyz VPSadenines Useaie ysiney aNd Uae adenine for each grid box in the interface grid WV Sysiney I5 4 Store the sum of the production and consumption rates of lysine and adenine respectively Lysine and adenine consumption rates are proportional to 9 where is lysine or K adenine concent
32. ution time at potential loss of accuracy Assuming fast dif fusion we may be able to adopt a larger grid spacing without significant loss of simulation accuracy We also try 20 and 40 um Interested readers may experiment with different grid spacings to find the optimal grid spacing void ModelRoutine updatelfGridSpacing REAL amp ifGridSpacing MODEL START ItGrid5pacing 5 0 7 5 0 20 0 or 40 0 to set the interface grid Spacing to 5 0 20 0 240 0 um MODEL END return Code 1 Model routine to set the interface grid spacing 2 updateOptModelRoutineCallInfo Set the number of rounds to update variables associated with the interface grid at the begin ning and at the end of a state and grid time step To set the model specific interface grid state variables properly Subheading 1 3 Bzocellton should be configured to invoke model routines updating interface grid state variables based on model specific rules once at the beginning of a state and grid time step and once at the end of the step in order to reset the sum of lysine consumption and production rates see Note 2 10 11 Simulating Microbial Community Patterning Using Biocellion 245 updateDomainBdryType Set domain boundary types Periodic boundary conditions are applied in the x and y directions A nonperiodic boundary condition is applied in the z direction and cells are not allowed to pass the upper and lower end of the
33. with a single 6 core micro processor Intel X5650 2 67 GHz Biocellion also accelerates larger higher resolution simulations using multiple compute nodes Biocellion runs on multicore PCs workstations cluster computers cloud computers and supercomputers This article pertains to mul ticore PCs Running on different systems does not require model code changes The current version of Bzocellion runs only on x86 compatible systems PCs with Intel or AMD microprocessors are x86 compatible A 64 bit Linux operating system needs to be 242 Seunghwa Kang et al O hours 100 hours 200 hours 300 hours 400 hours 500 hours 500 hours Fig 3 Yeast cell growth the top six figures viewed from an oblique angle from the top a 2D vertical cross section of the yeast colony the next figure and a 2D vertical cross section from the wet lab experiment the bottom most figure reproduced from Momeni et al s paper 5 Red cells are R strain cells green cells are G strain cells and black cells are dead cells installed on the target system Compiling Bzocellton model code requires the GNU gcc compiler pre installed in most Linux sys tems Intel icc compiler or some other C compiler we have tested only with gcc and 1cc When compiling model code users may set the check flag to verify their code or disable the check for 2 1 Installing Intel Thread Building Blocks 22 Installing Biocellion 3 Methods 3 1 Model Configuration
Download Pdf Manuals
Related Search