Home

Biocellion 1.1 User Manual

1. Interface Grid Spacing Interface Grid PDE Buffer Grid Spacing PDE Buffer Grid Figure 1 6 Simulation domain The simulation domain is decomposed into multiple partitions multiple compute processes on a cluster computer work on different sets of partitions Note that a single multithreaded compute process can exploit multiple cores on a compute node to work on a single partition Biocellion users set the partition size Biocellion outputs simulation results for different partitions in differ ent files If Biocellion restarts from checkpoint data Biocellion reads separate checkpoint files corresponding to distinct partitions Users can set each partition as a normal partition or a PDE buffer partition 1 3 Time Step Sizes Biocellion has three main computational modules with distinct time steps Communication inter vals also vary for different pairs of computational modules The baseline time step is the largest time step used by Biocellion Biocellion uses this time step size to simulate direct physico mechanical interactions between discrete agent pairs and to simulate discrete agent movement birth death and secretion The computational module simulating direct short range interactions communicates with the other two modules once every baseline time step Our assum
2. 4 5 22 updatePDEBufferDirichletBCVal 4 5 23 updatePDEBufferNeumannBCVal Simulation Output 4 6 1 updateSpAgentOutput 4 6 2 updateSummaryVar Biocellion Installation and Ezecution 6 1 6 2 6 3 6 4 6 5 6 6 6 7 6 8 System Requirements Compiling Biocellion 6 2 1 Compiling Model Code 6 2 2 Compiling the Biocellion Core Framework Simulation Configuration File 6 3 1 Required Elements 6 3 2 Optional Elements Execution on Desktop PCs and Workstations Execution on Clusters Execution on Amazon EC2 Debugging Support Visualization Using Paraview xii 53 53 53 54 54 54 55 55 55 57 59 59 59 59 59 60 60 61 62 62 63 63 63 6 8 1 Visualizing Discrete Agents 6 8 2 Visualizing Molecular Concentrations AP Tipsand aa AI A a ee A ele MR AE ee 8 0 C Basics Relevant to Biocellion 9 0 Bibliography xiii 63 64 65 67 69 1 1 1 2 1 3 1 4 1 5 1 6 1 7 4 1 4 2 4 3 4 4 List of Figures Discrete agent Shapes oi es oc ea e ea a eee et ede TE A junction between two discrete agents Uninhabitable region Setti
3. updatelfGridSpacing updateOptModelRoutineCallInfo 1X 31 32 32 32 32 32 32 33 33 33 33 33 34 34 34 34 34 35 37 37 38 38 38 4 3 4 2 3 updateDomainBdryType 4 24 updatePDEBufferBdryType 4 2 5 updateTimeStepinfo 4 2 6 updateSyneMethod 4 2 7 updateSpAgentinfo 4 2 8 updateJunctionEndinfo 4 29 updatePDEmfo 4 2 10 updatelfGridModelVarmnfo AD LL updateRNGINO yc na MDA E NESSE EEN 4 2 12 updateFileOutputinfo 4 2 13 updateSummaryOutputInfo 4 2 14 updateGlobalData ALIAS Ge Dolo Die A D eeh ep ASA AA A SE anga A E RE gl O Pd AE 42 IT SOLE DEUCE Ss RL DA AA RS A A Ee KU 4215 wena bia IA si is ii a e BOS AA Individual Agent Behavior 4 3 1 adSpAgenis 4 3 2 spAgeniCRNODERHS 4 3 3 updateSpAgentState 4 3 4 spAgentSecretionBySpAgent 4 3 5 updateSpAgentBirthDeath 4 3 6 adjustSpAgent 38 39 40 41 41 41 41 41 42 42 42 42 42 42 43 43
4. v_val Set the entire set of variables in the odeNetIdx ODE system to v_val v_val s size should coincide with the number of variables in the ODE system e REAL getModelReal const S32 idx const Return the value of the idx model specific REAL type variable e Vector lt REAL gt getModelRealArray void const Return a Vector array holding the entire set of model specific REAL type variable values void setModelReal const S32 idx const REAL val Set the value of the idx model specific REAL type variable to val void setModelRealArray const Vector lt REAL gt amp v val Set the entire set of model spe cific REAL type variables to v_val v val s size should coincide with the number of model specific REAL type variables for this discrete agent type e void incModelReal const S32 idx const REAL inc Increment the value of the idx model specific REAL type variable by inc e 32 getModellnt const S32 idx const Return the value of the idx model specific 32 type variable e Vector lt S32 gt getModelIntArray void const Return a Vector array holding the entire set of model specific S32 type variable values void setModellnt const S32 idx const S32 val Set the value of the idx model specific S32 type variable to val e void setModelIntArray const Vector lt S32 gt amp v_val Set the entire set of model specific S32 type variables to v_val v_val s size should coincide with the number of mode
5. 4928 gt sets the domain size to 128 x 128 x 4928 assuming the interface grid spacing Domain sizes in the x y and z directions should be equal to or larger than the minimum partition size which is 4 If adaptive mesh refinement AMR is used to solve PDEs domain sizes in the x y and z directions should be a positive integer multiple of AMR ratio AMR levels 1 x 4 See Section 6 3 2 for AMR ratio default value is 2 init_data Set the initialization method partition size and initialization file type and path Initialization method can be either code or file If the initialization method is set to file Biocellion users need to set file type binary is the only valid option at this point and path For example lt init data partition_size 32 src file file_type binar y path data input gt sets Biocellion to set partition size to 32 and initialize using check point data under data input Partition size should be equal to or larger than the minimum partition size which is 4 If adaptive mesh refinement AMR is used to solve PDEs partition size should be a positive integer multiple of AMR ratio AMR levels1 A See Section 6 3 2 for AMR ratio default value is 2 Note that the number of partitions in the simulation domain should be equal to or larger than the number of MPI processes in the target system to avoid idling MPI processes if the amount of computing in e
6. 43 43 43 44 44 45 45 4 4 4 5 4 3 7 divideSpAgent Physico Mechanical Interaction Between Agents 4 4 1 4 4 2 4 4 3 initJunctionSpAgent computeForceSpAgent computeExtraMechIntrctSpAgent State Changes in the Extra cellular Space 4 5 1 4 5 2 4 5 3 4 54 4 5 5 4 5 6 4 5 7 4 5 8 4 5 9 4 5 10 4 5 11 4 5 12 4 5 13 4 5 14 4 5 15 4 5 16 4 5 17 initlfGridVar updatelfGridVar updatelfGridKappa updatelfGridAlpha updatelfGridBetalnlfRegion updatelfGridBetaPDEBufferBdry updatelfGridBetaDomainBdry updatelfGridRHSLinear adjustIfGridRHSTimeDependentLinear updatelfGridRHSTimeDependentSplitting updatelfGridAMRTags updatelfGridDirichletBCVal updatelfGridNeumannBCVal initPDEBufferPhi updatePDEBufferKappa updatePDEBufferAlpha update PDEBufferBetalnPDEBufferRegion xi 45 46 46 46 47 47 47 48 48 49 49 49 50 50 50 51 51 51 52 52 52 52 53 5 0 6 0 4 6 Sample Implementations 4 5 18 updatePDEBufferBetaDomainBdry 4 5 19 updatePDEBufferRHSLinear 4 5 20 adjustPDEBufferRHSTimeDependentLinear 4 5 21 updatePDEBufferRHSTimeDependentSplitting
7. Termination 1 term M 4 2 5 updateTimeStepInfo void updateTimeStepInfo TimeStepInfo amp timeStepInfo Set baseline time step size and state and grid time step size This model routine updates timeStepInfo See Section 2 4 5 to find details about the TimeStepInfo class 40 4 2 6 updateSyncMethod void updateSyncMethod sync method e amp extraMechIntrctSyncMethod sync method e amp adjustG ridPhiAndModelVarSyncMethod Set the synchronization method when a single variable is updated in multiple model routine calls This model routine updates extraMechIntrctSyncMethod for extra direct mechanical interactions between agent pairs see computeExtraMechIntrctSpA gent in Section 4 4 3 and adjustGridPhi AndModelVarSyncMethod to update extra cellular space state variables based on model specific rules see updatelfGridVar in Section 4 5 2 See the explanation about sync_method_e in Sec tion 2 3 for the supported synchronization mechanisms 4 2 7 updateSpAgentInfo void updateSpA gentInfo Vector lt SpAgentInfo gt amp v_spAgentInfo Provide the information about the discrete agent types in the user model This model routine up dates v_spAgentInfo The size of v_spA gentInfo coincides with the number of discrete agent types Each vector element provides the information about a single discrete agent type See Section 2 4 7 for additional details about the SpA gentInfo class 4 2 8 update JunctionEndInfo voi
8. and z directions This is relevant when running Biocellion on a large cluster with a large number of nodes and hierarchical interconnection network Biocellion allows users to group closely located nodes Multiple partitions in a single super partition is processed by the MPI processes in a single node group and this reduces communication between distant nodes and improves the performance The sizes of a super partition in the x y and z directions should be a positive integer multiple of the partition size Note that the number of super partitions in the simulation domain should be equal to or larger than the number of node groups in the target system to avoid idling node groups if the amount of computing in each super partition varies significantly more super partitions are necessary to load balance lt super partition x 128 ya 128 za 128 gt sets the super partition sizes in the x y and z directions to 128 The default super 61 partition size in the x y or z direction is the smallest positive integer multiple of the partition size that is equal to or larger than the simulation domain size in the x y or z direction e interval Set the intervals for summary output load balancing AMR regridding and check point data output For example lt interval summary 1 load_balance 100 regriddi ng 100 checkpoint 600 gt sets Biocellion to print the summary every baseline time s
9. cons t UBAgentData amp ubAgentData const Vector lt REAL gt amp v_gridPhi const Vector lt REAL gt amp v gr idModelReal const Vector lt S32 gt 4 v_gridModellnt REAL amp gridBeta Set the diffusion coefficient for the elemidx molecular species for a unit box face between an interface grid partition and a PDE buffer partition The unit box face is orthogonal to a unit vector in the x if dim is 0 y if dim is 1 or z if dim is 2 direction vidx indexes the unit box in the interface grid partition ubAgentData holds the entire set of discrete agents in this unit box v_gridPhi v_gridModelReal and v_gridModelInt are the molecular concentrations model specific REAL type variables and model specific 32 type variables for the unit box in the interface grid partition respectively This model routine updates gridBeta 49 4 5 7 updatelfGridBetaDomainBdry void updatelfGridBetaDomainBdry const 32 elemlIdx const 32 dim const VIdx amp vIdx const U BAgentData amp ubAgentData const Vector lt REAL gt amp v_gridPhi const Vector lt REAL gt amp v_grid ModelReal const Vector lt S32 gt amp v_gridModellnt REAL amp gridBeta Set the diffusion coefficient for the elemIdx molecular species for a unit box face at the do main boundary The unit box face is orthogonal to a unit vector in the x if dim is 0 y if dim is 1 or z if dim is 2 direction vIdx indexes the unit box in the interface grid partition ubAgent Data hold
10. idx t type variables Public static member variables e VIdx ZERO A constant three dimensional zero vector i e 0 0 0 Vidz UNIT A constant three dimensional unit vector i e 1 1 1 e VIdxA BASIS 3 A set of three dimensional basis vectors A BASIS 0 1 0 0 A BASIS 1 0 1 0 and A_BASIS 2 0 0 1 Public non static member functions e x Arithmetic operators e l lt lt gt gt Comparison operators 15 e idx An access operator returns the idx element of the internal array holding three idx_t type variables idx should be 0 1 or 2 Public static member functions e S64 getIdx3Dto1lD const VIdx amp vide const VIdx amp regionSize Re turn vidx 0 x regionSize 1 x regionSize 2 vIdx 1 x regionSize 2 vIdx 2 e S64 getIdx3DtolD const VIdx amp a vIdx const idx_t size Returnvidx 0 x size x size vidx 1 x size vidx 2 e S64 getIdx3DTolD const idxt i const idxt j const idxt k const VIdx amp regionSize Return i x regionSize 1 x regionSize 2 j x regionSize 2 k e S64 getIdx3DTolD const idx_t i const idx_t j const idxt k const idx_t size Returni x size x size j x size k e VIdx getIdx1DTo3D const S64 idx const VIdx amp regionSize Per form the reverse operation of getIdx3DTo1D const Vidx amp vIdx const Vidx amp reg
11. state and the displacement from the original position disp state initially holds the state of this discrete agent disp should not ex ceed the interface grid spacing 4 3 7 divideSpAgent void divideSpA gent const VIdx amp vIdx const AgentJunctionInfo amp junctionInfo const VReal amp vO Deet const AgentMechIntrctData amp mechIntrctData const Vector lt NbrBox lt REAL gt gt amp v gridP hiNbrBox const Vector lt NbrBox lt REAL gt gt amp v_gridModelRealNbrBox const Vector lt NbrBox lt 32 gt gt amp v_gridModelIntNbrBox SpAgentState amp motherState VReal amp motherDisp SpAgentS tate amp daughterState VReal amp daughterDisp Vector lt BOOL gt v_junctionDivide BOOL amp moth erDaughterLinked JunctionEnd amp motherEnd JunctionEnd amp daughterEnd 45 Divide this discrete agent to a mother discrete agent and a daughter discrete agent vidr indexes the unit box this discrete agent is located in junctionInfo stores the list of junctions between this discrete agent and its neighboring discrete agents vOffset is the offset from the center of the unit box this discrete agent belongs to mechIntrctData stores the temporary mechanical interaction state of this discrete agent v_gridPhiNbrBox v_gridModelRealNbrBox and v_gridModelIntNbr Box provide the information molecular concentrations model specific REAL type attributes and model specific S32 type attributes respectively about the extra cellular space for this
12. 24 135 MG Model Vari AA o RER EE EN 2 4 14 RNGInfo 2 4 15 FileOutputInfo 2 4 16 SummaryOutputlnfo 2 4 17 SpAgentState 2 4 18 JunctionEnd 2 4 19 AgentJunctionInfo 2 4 20 SpAgent 2 4 21 UBAgentData 2 4 22 ErtraMechintretData 2 4 23 AgentMechlntrctData 2 4 24 IfGridUpdate 2 4 25 NbrBox 2 4 26 IfGridBoxData 3 0 Biocellion Support Functions 3 1 Simulation Instance Information 3 1 1 getCurBaselineTimeStep 3 1 2 getCurStateAndGridTimeStep 3 1 3 getCurPDETimeStep 3 1 4 getGlobalDataRef 20 20 21 22 22 22 23 23 24 25 26 27 27 28 28 28 29 31 31 31 31 31 31 4 0 3 1 10 3 1 11 3 1 12 3 2 Random Number Generator 32 1 3 3 Extra Support Functions 3 3 1 3 3 2 3 3 3 Biocellion Model Routines 4 1 Model Routine Invocation Timings 4 2 Model Configuration 4 2 1 42 2 getRecentSummaryRealVal getRecentSummarylntVal getDomainSize getPartSize getSimlnitlype getParticleNumExtraOutputVars getModelParam getAMRRatio getModelRand computeSphereUB VolOvlpRatio changePosFormat2LvTolLv changePosFormatlLvTo2Lv
13. Click the Display tab Set Color by to color Click Zoom To Data This maps discrete agents to a set of spheres Using the Mask Points filter instead of the Glyph filter maps a discrete agent to a point instead of a sphere this is a computationally less expensive option to visualize discrete agent data Visualizing Molecular Concentrations Click the File menu and select the Open item Open a set of files in the VTK hierarchical box data files format vtm When a pop up menu appears select the VTK Hierarchical Box Data Files format Click the Properties tab at the bottom left side and click the Apply icon Select the Filters menu the Alphabetical sub menu and the Slice item Select the Property tab and adjust the plane to slice and click the Apply icon Select the Display tab Set Color by to phi Click Zoom To Data This displays the molecular concentrations at the slice plane 64 7 0 Tips and Caveats We strongly encourage Biocellion users to enable the checks on model routine outputs and input arguments to Biocellion support routines Section 6 2 1 if the simulation does not work properly This often guides users to find errors on their model routines The relative and absolute error norm thresholds for multigrid iterations Section 2 4 10 have significant impact on the accuracy and speed of simulation If the
14. National Technical Information Service U S Department of Commerce 5285 Port Royal Rd Springfield VA 22161 ph 800 553 6847 fax 703 605 6900 email ordersOntis fedworld gov online ordering http www ntis gov ordering htm This document was printed on recycled paper 9 2003 Acknowledgments Support for this research was provided by the Extreme Scale Computing Initiative and the Funda mental 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 ill Changes Changes From Version 1 0 e Modify the Biocellion API related to file output Added the update FileOutputInfo Section 4 2 12 model routine and the FileOutputInfo class Section 2 4 15 Removed the fileOutput member variable from the GridPhilnfo class Section 2 4 9 updateFileOutputInfo controls file output for molecular concentrations in the extracel lular space Removed num extra from the output element Section 6 3 1 updateFileOutputinfo controls file output for particles Added start_x start_y start_z size_x size_y and size_z to the output element Sec tion 6 3 1 Users can generate file output data for only a sub region of the simulation domain e Modify the Biocellion API related to sum
15. Set the number of sockets to the number of sockets per shared memory node for non uniform memory access NUMA binding Set the number of sockets per node to 1 to disable NUMA binding For example lt system num_node_groups 1 num_nodes_per_group l num socket s_pernode 1 max load_imbalance 1 05 num_threads 8 gt sets the number of node groups the number of nodes per group and the number of sockets per node to 1 sets Biocellion to perform load balancing in the designated interval set in the system configu ration file see the explanation about the interval element in this section if the MPI process with the highest load has more than 1 05 times load over the average and sets the num ber of threads per MPI process to 8 The default number of node groups is 1 The default value for the number of nodes per node group is equal to the number of MPI processes in the MPI communicator The default value for the number of sockets per node is 1 disable NUMA binding The default threshold for load balancing is 1 2 and the default number of threads per MPI process is the number of cores or hardware threads in a system with hard ware multi threading support in the target system Note that the number of node groups x the number of nodes per group x the number of sockets per node should coincide with the number of MPI processes in the MPI communicator e super partition Set the super partition size in the x y
16. a model specific REAL variable IF_GRID_VAR_TYPE_MODEL_INT Set the type to a model specific 32 variable e sim init type e Used to set the simulation initialization method SIM INIT TYPE CODE Initialize inside the model routines SIM_INIT_TYPE_BINARY Initialize using checkpoint data e particle output type e Used to set the output file format for discrete agents PARTICLE_OUTPUT_TYPE_PTVU Set the output file format to the Paraview VTK partitioned unstructured grid file format ptvu This is the only valid option at this point e grid output type e Used to set the output file format for grid values GRID_OUTPUT_TYPE_VIM Set the output file format to the VTK hierarchical box data files format vtm This is the only valid option at this point summary type e Biocellion supports a summary reduction mechanism for interface grid variables Users add a fixed number of grid state variables to every unit box in the interface grid or assign a fixed number of attributes to the interface grid and for each attribute users can ask Biocellion to visit every unit box to reduce the variables for the attribute to a single value summary type e is used to set the reduction method SUMMARY TYPE SUM Set to print the sum of the interface grid variables SUMMARY TYPE AVG Set to print the average of the interface grid variables SUMMARY
17. face e ode stiff_e Inform the stiffness of the ODE system to solve Biocellion uses different mathematical algorithms based on the provided information using the Intel ODE solver Intel Corporation 2008 ODE_STIFF_NORMAL Indicate that the ODE system has low or medium stiffness The dodesol_rkm9st function in the Intel ODE library which is based on the fourth order Merson s method and the first order multistage method of up to and in cluding 9 stages with stability control Intel Corporation 2008 is used to solve the ODE system ODE STIFF HIGH Indicate that the ODE system is stiff The dodesol mk521fn function in the Intel ODE library which uses the implicit method based on L stable 5 2 method with the numerical Jacobi matrix which is computed by the routine Intel Cor poration 2008 is used to solve the ODE system e sync method e Biocellion users can update temporary mechanical interaction state vari ables of a discrete agent in multiple function calls evaluating direct interactions between the discrete agent and the neighboring discrete agents one function call per pair of discrete agents Biocellion provides a model routine to edit grid state variables using model specific rules and the routine is called once per unit box in the interface grid A single invocation of the model routine can update variables associated with an interface grid unit box and its neighboring 26 unit boxes A
18. norm thresholds are too large users may see no change at all in molecular concentrations even though the source terms in PDEs are set to non zero values Setting these values to unnecessarily small values can significantly slow down the simulation on the other hand Be cautious when setting these parameters Note that double instead of REAL is used in spA gentCRNODERHS this is mainly be cause the Intel ODE solver does not support single precision arithmetic Be cautious about setting Neumann boundary values when k see Sections 4 5 3 and 4 5 15 is set to a value smaller than 1 0 Biocellion does not consider x in setting boundary values x is used only in computing diffusion flux For example users may expect there is no net diffusion if the Neumann boundary value for a boundary unit box face is set to 0 0 Biocellion sets the molecular concentration of the box outside the simulation domain to the molecular concentration of the box inside the simulation domain if the Neumann boundary value is 0 0 the gradient ignoring K becomes 0 0 Biocellion assumes that k outside the simulation domain is 1 0 and if the k value for the box inside the simulation value is smaller than 1 0 JI can have a non zero value Users may set the diffusion coefficient at the boundary to 0 0 instead to apply the zero flux boundary condition 65 8 0 C Basics Relevant to Biocellion We provide several links for users who are new to C e http www c
19. than the state and grid time step by partitioning a single state and grid time step to a positive integer number of PDE time steps To solve the splitting PDE Biocellion solves ODEs for the reaction part for the duration of a single PDE time step first and sets the PDE source term based on the changes Then the resulting linear reaction diffusion PDE is solved for one or more diffusion time steps a single PDE time step can be further partitioned to multiple diffusion time steps when the splitting scheme is used Biocellion uses a variable time step size based on error estimation to solve ODEs following the scheme used in the Intel ODE solver and Biocellion users set the minimum ODE time step size The single PDE time step size becomes the maximum time step size to solve ODEs 1 5 Future Roadmap Over time we plan to update Biocellion to accommodate an ever broadening spectrum of biological system models In the short term we plan to extend Biocellion to support mapping a single cell to multiple discrete agents and to provide advection reaction diffusion PDE solvers Mapping a cell to more than one discrete agent will allow modeling of cell shapes and mechanical interactions between cells with higher fidelity For example mapping a budding yeast cell to two discrete agents more realistically reproduces the shape of a real budding yeast cell Figure 1 7 Figure 1 7 Modeling budding yeast using two discrete agents per cell The budding yea
20. this discrete agent due to direct mechanical interactions with neighboring discrete agents and the entire set of extra temporary mechanical in teraction state variables Public non static member variables e VReal force The sum of the forces on this discrete agent e ExtraMechIntrctData extraMechIntrctData This variable captures the model specific extra temporary mechanical interaction state 2 4 24 IfGridUpdate This class instance holds the information required to update a single state variable of a unit box in the interface grid e if_grid_var_type_e type This variable sets the type of the interface grid state variable to update e S32 elemldx This variable sets the index of the interface grid state variable to update union REAL realVal S32 intVal val This variable stores the new value of the interface grid state variable to update 2 4 25 NbrBox This template class instance holds a set of 27 variables of type T type T is given when this class instance is defined for an interface grid unit box and its 26 neighboring boxes one value per box These variables are valid only for the boxes in a valid interface region An offset variable in each dimension xOffset yOffset zOffset vOffset 0 vOffset 1 or vOffset 2 should be 1 0 or 1 xOffset yOffset zOffset or alternatively vOffset 0 vOffset 1 vOffset 2 0 O 0 indexes the center box e BOOL getValidFlag const S32 xOffset const S32 yOffs
21. x 1 0 to turn this check off 19 e BOOL setNegToZero If setNegToZero is set to true Biocellion automatically resets nega tive PDE variables equal to or larger than errorThresholdVal to 0 0 Set to false to turn this feature off 2 4 10 MGSolvelInfo This class instance holds the model specifics about the PDE multigrid solver Public non static member variables e S32 numPre Set the number of smoothing steps before averaging in the multigrid solver Set to 3 unless you wish to tune the multgrid solver based on your specific PDE problems e S32 numPre Set the number of smoothing steps after averaging in the multigrid solver Set to 3 unless you wish to tune the multgrid solver based on your specific PDE problems e S32 numBottom Set the number of smoothing steps at the bottom level Set to 3 unless you wish to tune the multgrid solver based on your specific PDE problems e BOOL vCycle Set to true to use V cycle or set to false to use W cycle in the multigrid solver Set to true unless you wish to tune the multgrid solver based on your specific PDE problems e S32 maxlters Set the maximum number of V cycles or W cycles Se to 30 unless you wish to tune the multigrid solver based on your specific PDE problems e g if your PDE converges very slowly a larger value could be more appropriate If the PDE multigrid solver does not reach convergence within maxlters iterations the solver exits prematurely and prints a warni
22. 2012 Intel Corporation 2008 for additional details The steady state PDE and the time dependent linear PDE update concentrations of a single molec ular species The time dependent PDE based on the splitting scheme update the concentrations of one or more molecular species reacting in the extra cellular space Biocellion supports adaptive mesh refinement AMR to solve PDEs using CHOMBO Colella etal 2012 The grid spacing used for the interface grid is the finest spacing in the AMR hierarchy Users can set the number of AMR levels and select either 2 or 4 as the refinement ratio between two consecutive levels Users mark the desired AMR level for each unit box in the interface grid The coarsest AMR level is assumed for the PDE buffer The actual AMR grid may use a finer grid spacing than the grid spacing specfied by the user for some regions to satisfy necessary conditions to maintain numerical efficiency e g to limit the number of boxes because processing a large number of tiny boxes can be computationally inefficient and correctness e g see the proper nesting condition in the CHOMBO manual A variant of the second order accurate Lo stable Runge Kutta scheme see the CHOMBO manual and the Twizell et al s paper Twizell et al 1996 for additional details is used for time stepping This often allows Biocellion users to adopt a larger time step size than explicit methods for stiff PDEs A single PDE time step can be set to a smaller value
23. Biocellion redefines the following basic data types to promote portability across different plat forms and to support both single precision and double precision floating point arithmetic without code modification e REAL Either single or double precision floating point number float or double based on compile time setup e U8 U16 U32 U64 8 16 32 and 64 bit unsigned integers respectively e S8 S16 S32 S64 8 16 32 and 64 bit signed integers respectively BOOL Boolean data type true or false BOOL replaces the C bool data type Biocellion uses BOOL instead of bool to prevent the C standard template library STL vector from packing multiple boolean variables up to 8 into a single byte Such packing reduces memory consumption but requires additional bitwise operations to access and modify vector elements If vector elements are concurrently updated by multiple threads updating two different elements packed into the same byte requires locking or using other data synchronization primitives we expect users to write only sequential code so data synchronization is not an issue for Biocellion users but we define the BOOL type to avoid additional data type conversion overhead between user routines and the Biocellion framework code 2 2 Extra Scalar Data Types Biocellion defines additional scalar data types to maintain portability across different platforms and save memory e agentType t 16 bit unsigned integer by de
24. CELLION_ROOT framework main biocellion DP SPAGENT D EBUG is a Biocellion executable compiled with debugging support and biocellion DP SPA GENT OPT is an executable that delivers high performance DP stands for double precision use the executables with SP instead of DP if single precision floating point arithmetic is sufficient Users execute Biocellion by typing biocellion DP SPAGENT OPT simulation configuration file 6 5 Execution on Clusters Use biocellion DP SPAGENT MPI OPT or biocellion DP SPAGENT MPI DEBU G for debugging support to run on clusters with the Ethernet interconnect and MPICH2 http www mpich org For clusters with the InfiniBand Interconnect or other proprietary in terconnects or different flavors of MPI libraries users need to recompile the Biocellion framework using Global Arrays and MPI libraries configured for the target system mpiexec n numb er of MPI processes machinefile machine file name biocellion DP SPAGNT MPI OPT simulation configuration file name launches multiple MPI processes to run a large simulation using multiple nodes Refer to MPI manuals for additional details To run Biocellion on a cluster with a job scheduler refer to a job scheduler manual or contact a system administrator 62 6 6 Execution on Amazon EC2 Launch EC2 instances we use either CC1 Cluster Compute cc1 4xlarge or CC2 Cluster Com pute
25. IFFUSION TIME DEPENDENT LINEAR Use time stepping to find a solution of a reaction diffusion linear PDE PDE TYPE REACTION DIFFUSION TIME DEPENDENT SPLITTING Use time stepping to find a solution of a reaction diffusion PDE The time splitting technique Strang 1968 is applied to simulate reactions among multiple molecular species in the extra cellular space 12 PDE_TYPE_ADVECTION_REACTION_DIFFUSION_TIME_DEPENDENT_LINEAR Re served for future use PDE TYPE ADVECTION REACTION DIFFUSION TIME DEPENDENT SPLITTING Reserved for future use e bc type e Used to set PDE boundary condition types BC TYPE DIRICHLET CONST Apply Dirichlet boundary condition with a constant boundary value for the entire region of the domain boundary face at the low or high side in the x y or z direction BC_TYPE_ DIRICHLET_MODEL Apply Dirichlet boundary condition and boundary values are set in the granularity of a boundary grid unit box face BC TYPE NEUMANN CONST Apply Neumann boundary condition with a constant boundary value for the entire region of the domain boundary face at the low or high side in the x y or z direction BC TYPE NEUMANN MODEL Apply Neumann boundary condition and boundary val ues are set in the granularity of a boundary grid unit box
26. PNNL XXXXX Biocellion 1 1 User Manual Seunghwa Kang and Simon Kahan March 31 2014 Prepared for the U S Department of Energy under Contract DE AC05 76RL01830 DISCLAIMER United States Government Neither the United States Government nor any agency thereof nor Battelle Memorial Institute nor any of their employees makes any warranty express or implied or assumes any legal liability or responsi bility for the accuracy completeness or usefulness of any information ap paratus product or process disclosed or represents that its use would not infringe privately owned rights Reference herein to any specific commercial product process or service by trade name trademark manufacturer or other wise does not necessarily constitute or imply its endorsement recommendation or favoring by the United States Government or any agency thereof or Battelle Memorial Institute The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof PACIFIC NORTHWEST NATIONAL LABORATORY operated by BATTELLE for the UNITED STATES DEPARTMENT OF ENERGY under Contract DE ACO5 76RLO1830 Printed in the United States of America Available to DOE and DOE contractors from the Office of Scientific and Technical Information P O Box 62 Oak Ridge TN 37831 0062 ph 865 576 8401 fax 865 576 5728 email reports adonis osti gov Available to the public from the
27. RHSLinear PB 16 adjustPDEBufferRHSTimeDependentLi near PB 17 updatePDEBufferRHSTimeDependent Splitting PB and grid time steps 2 computeExtraMechIntrciSpAgent AP Yes i L 1 updatelfGridVar IB PR Simulate mechanical interactions Update grid values pre E e SS ES e Update agent states Solve PDEs Update grid values post 1 updatelfGridVar 1B 18 updatePDEBufferDirichletBCVal PF 19 updatePDEBufferNeumannBCVal PF Nae Simulate agent birth and death 1 updateSpAgentOutput A File 8 summary output 2 updateSummaryVar IB File gisummary output amp summary File gisummary output A ewech E 1 addSpAgents P 1 computeForceSpAgent AP 1 spAgentCRNODERHS A 2 updateSpAgentState A spAgentSecretionBySpAgent A updateSpAgentBirthDeath A adjustSpAgent A divideSpAgent A disappearSpAgent A Sai 1 updatelfGridAMRTags IB End main loop Figure 4 2 Main loop 4 2 4 updatePDEBufferBdryType void update PDEBufferBdryTypel pde buffer bdry type e amp pdeBufferBdryType Set the boundary type between an interface grid partition and a PDE buffer partition This model routine updates pdeBufferBdryType See the explanation about pde_buffer_bdry_type_e in Sec tion 2 3 for the supported PDE buffer boundary types 39 updatelfGridKappa IB updatelfGridAlpha IB updatelfGridBetaln
28. TYPE MIN Set to print the minimum value of the interface grid variables 14 SUMMARY TYPE MAX Set to print the maximum value of the interface grid variables e rng_type_e Biocellion provides random number generators internally using Intel math kernel library MKL mg type e is used to set the random number generator type RNG TYPE UNIFORM Set the random number generator to generate random numbers with the uniform distribution BNG TYPE GAUSSIAN Set the random number generator to generate random num bers with the Gaussian distribution BNG TYPE EXPONENTIAL Set the random number generator to generate random numbers with the exponential distribution BNG TYPE GAMMA Set the random number generator to generate random numbers with the gamma distribution 2 4 Classes 2 4 1 Vector Vector is a wrapper class of the C standard template library STL vector STL vector is a data type for a variable size array We create this wrapper to perform array index bound checks sim ilar to CHOMBO Colella et al 2012 without re compiling both the Biocellion core framework and a model library Users can enable the boundary check for a model library without enabling the check for the core framework Visithttp www cplusplus com reference vector vector for additional details about STL vector 2 4 2 Vidz Vidz is used to index a grid box in the simulation domain VIdx has an internal array of three
29. The ODE solver cannot ensure required accuracy if epsilon is smaller than le 9 e REAL threshold If the absolute value of the solution is larger than threshold relative error is controlled Otherwise absolute error threshold x epsilon is controlled e REAL errorThresholdVal If any of the ODE system variables becomes smaller than er rorThresholdVal Biocellion prints an error message and aborts the program execution er rorThresholdVal should be equal to or smaller than 0 0 Set errorThresholdVal to REAL MAX defined in SBIOCELLION ROOT include base type h x 1 0 to turn this check off e REAL warningThresholdVal If any of the ODE system variables becomes smaller than warningThresholdVal but still no variable is smaller than errorThresholdVal Biocellion prints a warning message but does not abort the program execution warningThresholdVal should be equal to or smaller than 0 0 and should be equal to or larger than errorThreshold Val Set warningThresholdVal to REAL MAX defined in SBIOCELLION ROOT inclu de base type h x 1 0 to turn this check off e BOOL setNegToZero If setNegToZero is set to true Biocellion automatically resets nega tive ODE variables equal to or larger than errorThresholdVal to 0 0 Set to false to turn this feature off 2 4 5 TimeStepInfo This class instance holds the information required to set the sizes of baseline time step and state and grid time step Public non
30. X Wang W Liu Y Hao G Cai and YW Han 2013 Fusobacterium nucleatum Promotes Colorectal Carcinogenesis by Modulating E Cadherin Catenin Signaling via its FadA Adhesin Cell Host amp Microbe Sandersius SA CJ Weijer and TJ Newman 2011 Emergent cell and tissue dynamics from subcellular modeling of active biomechanical processes 69 Singh JS P Abhilash H Singh RP Singh and D Singh 2011 Genetically engineered bacteria An emerging tool for environmental remediation and future research perspectives Gene 480 1 2 1 9 Stamatakos GS 2010 Multiscale Cancer Modeling Ch In Silico Oncology Part I Clinically Oriented Cancer Multilevel Modeling Based on Discrete Event Simulation CRC Press Strang G 1968 On the construction and comparison of difference schemes SIAM Journal on Numerical Analysis 5 3 506 517 Twizell EH AB Gumel and MA Arigu 1996 Second order Lo stable methods for the heat equation with time dependent boundary conditions Advances in Computational Mathematics 6 333 352 Weinberg RA 2007 The biology of cancer Garland Science Xavier JB C Picioreanu and MCM van Loosdrecht 2005 A framework for multidimen sional modelling of activity and structure of multispecies biofilms Environmental Microbiology 7 8 1085 1103 70
31. ach partition varies significantly more partitions are necessary to load balance output Set the output path interval and file formats pvtu is the only valid option for discrete agent and vtm is the only valid option for molecular concentrations at this point For example lt output path data output interval 100 particle pvtu grid vim gt sets the file output path to data output the file output interval to once per 60 every 100 baseline time steps the particle output format to the parallel VTK unstructured grid format and the molecular concentration output format to the VTK hierarchical box data format Setting interval to 0 turns off file output 6 3 2 Optional Elements e model Set the model specific parameter string Model routines can access the parameter string using getModelParam Section 3 1 11 For example lt model param test gt sets the model specific parameter to test The model specific parameter is set to by default e verbosity Set the standard output verbosity 0 to 5 5 for the highest level of verbosity For example lt stdout verbosity 3 gt sets the standard output verbosity level to 3 The default verbosity level is 1 e system Set the number of node groups number of nodes per node group number of sock ets per node threshold to perform load balance and the number of threads per MPI pro cess
32. annBCVal const S32 elemidx const VReal amp pos const 32 dim const B OOL lowSide const REAL a_gridPhi 3 const Vector lt Vector lt REAL gt gt amp vv_gridModelReal const Vector lt Vector lt 8S32 gt gt amp vv_gridModellnt REAL amp bcVal Set the model specific Neumann boundary value for the elemIdx molecular species for a unit box face at the domain boundary This model routine is valid only when BC TYPE NEUMANN MODEL is selected pos locates the center of the unit box face The unit box face is orthogonal to a unit vector in the x if dim is 0 y if dim is 1 or z if dim is 2 direction The face is at the low side if lowSide is true or the high side if lowSide is false of the simulation domain a gridPhi holds molecular concentrations for the elemIdx molecular species for the three interface grid unit boxes closest to the boundary unit box face along the direction orthogonal to the face vv_grid ModelReal and vv gridModelInt hold model specific REAL and S32 type variables for the three boxes closest to the boundary unit box face This model routine updates bcVal 4 5 14 nitPDEBufferPhi void initPDEBufferPhi const Vidx amp startVIdx const Vidx amp pdeBufferBoxSize Vector lt REAL gt amp v_gridPhi Initialize molecular concentrations of a PDE buffer unit box startVIdx points the low end cor ner of the PDE buffer unit box assuming the interface grid spacing assuming that the interface grid is applied
33. artition and regionSize is the size of this partition ifGridHabitableBoxData in forms whether a unit box in this partition is habitable or not It is an error to add a discrete agent to an uninhabitable unit box This model routine updates v_spAgentVIdx v_spAgentState and v_spA gentOffset The sizes of v spAgentVIdx v_spAgentState and v_spAgentOffset should be same one vector element for one newly added discrete agent An element of v_spA gentVIdx lo cates the unit box to place a discrete agent An element of v_spAgentState stores the state of a discrete agent An element of v_spA gentOffset sets the discrete agent position offset within a unit box from the center of the unit box 4 3 2 spAgentCRNODERHS void spAgentCRNODERHS const S32 odeNetldx const Vidx amp vIdx const SpAgent amp spAgent co nst Vector lt NbrBox lt REAL gt gt amp v_gridPhiNbrBox const Vector lt NbrBox lt REAL gt gt amp v_grid ModelRealNbrBox const Vector lt NbrBox lt 832 gt gt amp v_gridModelIntNbrBox const Vector lt dou ble gt amp vy Vector lt double gt amp vf Set the right hand side the derivatives of the odeNetldx ODE system for this discrete agent vidx indexes the unit box this discrete agent is located in spAgent provides the information about 43 this discrete agent v_gridPhiNbrBox v_gridModelRealNbrBox and v gridModelIntNbrBox pro vide the information molecular concentrations model specific REAL type attributes and model sp
34. ary OutputInfo and updateSummaryVar model routines to Figures 4 1 and 4 2 Removed the getMGEllipticNormThreshold and getMGParabolicNormThreshold functions Section 3 1 Added the getSimInitType function Section 3 1 9 Removed the disappearSpAgent model routine Section 4 3 and the IfGridUpdate class Section 2 4 Added getRecentSummaryRealVal and getRecentSummaryIntVal functions Sections 3 1 5 and 3 1 6 respectively to access recent summary output values Replaced minVal in the ODENetInfo and GridPhilnfo classes Sections 2 4 4 and 2 4 9 re spectively with errorThresholdVal and warningThresholdVal vi Acknowledgments Contents Changes IA BE Eege AE ee ee e EA A 1 0 Introduction 1 1 Simulation Components 1 2 Simulation Domain 1 3 Time Step Sizes 1 4 Partial Differential Equation PDE Solvers 1 5 Future Roadmap 2 0 Biocellion Data Types 2 1 Basic Data Types 22 ERAS calar Data Eypes ts Des SRS See Se E BS OR e BS Sad 2 3 Enumerations 2 4 Classes 2 4 1 2 4 2 2 4 3 2 4 4 2 4 5 2 4 6 2 4 7 2 4 8 2 4 9 TimeStepInfo OptModelRoutineCalllnfo SpAgentInfo JunctionEndInfo 0 0 0 00 A GridPhilnfo vii 111 11 11 11 12 15 15 15 16 17 17 18 18 19 19 2 4 10 MGSolvelnfo 2 4 11 SplittingInfo 24 12 PDEInfo
35. box and the 26 neighboring boxes This model routine updates the states of the mother discrete agent and the daughter discrete agent motherState and daughterState respectively the displacements from the original position for the mother and daughter discrete agents motherDisp and daughterDisp re spectively and a junction between the mother and daughter discrete agents motherState initially stores the state of this discrete agent Users may overwrite the initial state v_junctionDivide in dicates a subset of the junctions if the v_junctionDivide element for a junction is set to true the mother agent retains the junction the mother agent will retain the daughter agent inherits the remaining junctions A junction is formed between the mother and daughter discrete agents if motherDaughterLinked is set to true with motherEnd and daughterEnd as the junction ends for the mother and daughter sides respectively motherEnd and daughterEnd are irrelevant if moth erDaughterLinked is set to false 4 4 Physico Mechanical Interaction Between Agents 4 4 1 initJunctionSpAgent void initJunctionSpA gent const Vidx amp vIdx0 const SpAgent amp spAgent0 const VIdx amp vIdx1 con st SpAgent amp spAgentl const VReal amp dir const REAL amp dist BOOL amp link JunctionEnd amp endo JunctionEnd amp end Initialize a junction between a discrete agent pair within the maximum direct interaction distance vidx0 indexes the unit box spAgent0 is locat
36. cc2 8xlarge using the Biocellion machine image refer to the Biocellion homepage http biocellion org for the latest machine image name Refer the Amazon EC2 homepage nttp aws amazon com ec2 for instructions to launch EC2 instances Log in to one of the launched instances Biocellion is pre installed under home ec2 user To run ona single EC2 instance refer to the instructions in Section 6 4 Intel TBB is pre installed and the LD_LIBRARY_PATH Linux environment variable is already set To run on multiple instances update mpd hosts file under home ec2 user to list the launched EC2 instance addresses Type mpdboot n number_of_EC2_instances f mpd hosts to launch an MPI daemon Type mpiexec n number of MPI processe s biocellion DP SPAGNT MPI OPT simulation_configuration_file_nam e to execute Biocellion using multiple instances Note that each instance has its own home e c2 user directory Updating a file under home ec2 user in one instance will not update the file with the same name under home ec2 user of different instances Create an Amazon S3 storage to share files across multiple EC2 instances 6 7 Debugging Support Biocellion installation includes binaries with debugging support Users can enable checks on model routine outputs and input arguments to Biocellion support functions by enabling the check flag on compile time as stated in Section 6 2 1 Users can set Biocellion to monitor ordinary differentia
37. cell movement paths Resat et al 2011 Biocellion provides two different mechanisms to update environment state variables e g molecu lar concentrations or other extra cellular space attributes One is editing state variables by model specific rules and the other is by solving partial differential equations PDEs Biocellion users also mark sub regions of the simulation domain as uninhabitable and discrete agents cannot pen etrate into an uninhabitable region Figure 1 3 Displacement set inside the model routine Inhabitable Region Actual displacement Figure 1 3 A discrete agent cannot move into an uninhabitable region Biocellion considers only the position of a discrete agent in computing the displacement and does not internally use other discrete agent attributes e g radius By default Biocellion al lows any discrete agent of non zero volume to overlap uninhabitable regions If this is undesirable the modeler can include code to test the agent s physical position and orientation adjusting the displacement so as to move the agent fully out of uninhab itable regions Biocellion has three main computational modules to simulate the above simulation components The three modules 1 update the state of individual agents 2 simulate short range physico mechanical interactions and 3 update the state of the extracellular environment Computa tional modules communicate with other modules to simulate biological syst
38. core framework to a different flavor of MPI library 59 Compiling the core framework is more involved and requires the Intel Thread Building Block library version 4 2 or later Intel math kernel library MKL Intel ODE library PNNL Global Arrays http www emsl pnl gov docs global and MPI in addition Compiling the core framework also requires a C compiler that supports c 0x extension recent versions of GNU gcc and Intel icc support the extension Users set DIST MEM PAR in SBIOCELLION_ROOT Makefile framework to yes to cre ate a Biocellion executable that runs using multiple MPI processes or to no to create a Biocellion version that runs using only one multi threaded process Users also need to update other direc tory paths in the file Typing make under BIOCELLION_ROOT framework places a newly created executable under SBIOCELLION ROOT framework main 6 3 Simulation Configuration File Biocellion asks users to provide specifics of a simulation instance in an XML file sample files are located under SB OCELL ON ROOT framework main 6 3 1 Required Elements e time step Set the number of baseline time steps to execute For example lt time step num baseline 60000 gt asks Biocellion to run for 60 000 baseline time steps domain Set the domain sizes in the x y and z directions For example lt domain x 128 y 128 z
39. d after the ODE systems for the discrete agent are updated 4 3 4 spAgentSecretionBySpAgent void spAgentSecretionBySpAgent const VIdx amp vIdx const AgentJunctionInfo amp junctionInfo con st VReal amp vOffset const AgentMechIntrctData amp mechIntrctData const Vector lt NbrBox lt REAL gt gt amp v_gridPhiNbrBox const Vector lt NbrBox lt REAL gt gt amp v gridModelRealNbrBox const Ve ctor lt NbrBox lt S32 gt gt amp v_gridModelIntNbrBox SpAgentState amp state Vector lt SpAgentState gt amp v_spAgentState Vector lt VReal gt amp v_spAgentDisp Secrete new discrete agents from this discrete agent v dx indexes the unit box this discrete agent is located in junctionInfo stores the list of junctions between this discrete agent and its neigh boring discrete agents vOffset is the offset from the center of the unit box this discrete agent be longs to mechIntrctData stores the temporary mechanical interaction state of this discrete agent v_gridPhiNbrBox v_gridModelRealNbrBox and v_gridModelIntNbrBox provide the information molecular concentrations model specific REAL type attributes and model specific S32 type at tributes respectively about the extra cellular space for this box and the 26 neighboring boxes This model routine updates this discrete agent s state state and the states v_spAgentState and the displacements from this discrete agent v_spAgentDisp of newly secreted discrete agents The displacement valu
40. d updateJunctionEndInfo Vector lt JunctionEndInfo gt amp v_junctionEndInfo Provide the information about the junction end types in the user model This model routine up dates v_junctionEndInfo The size of v_junctionEndInfo coincides with the number of junction end types Each vector element provides the information about a single junction end type See Section 2 4 8 for additional details about the JunctionEndInfo class 4 2 9 updatePDEInfo void updatePDEInfo Vector lt PDEInfo gt amp v pdelnfo Provide the information about the molecular species in the user model This model routine up dates v_pdelnfo The size of v_pdelnfo coincides with the number of molecular species Each vector element provides the information about a single molecular species See Section 2 4 12 for additional details about the PDEInfo class 4 2 10 updatelfGridModelVarInfo void updatelfGridModelVarInfo Vector lt IfGridModelVarInfo gt amp v_ifGridModelRealInfo Vecto r lt IfGridModelVarInfo gt amp v_ifGridModelIntInfo Provide the information about the model specific attributes for the extra cellular space in the user model This model routine updates v_ifGridModelReallnfo for REAL type attributes and v_if 4 GridModelIntInfo for S32 type variables The sizes of v ifGridModelRealInfo and v_ifGridMod elIntinfo coincide with the number of REAL type and S32 type model specific attributes respec tively Each vector element provides the informatio
41. dateSpAgentOutput const Vidx amp vidx const SpAgent amp spAgent REAL amp color Vector lt REAL gt v_extra Set the output variables for a discrete agent vidr indexes the unit box this discrete agent spA gent belongs to This model routine update the color variable color and the model specific extra output variables v_extra The number of model specific output variables for a discrete agent is set inside the simulation configuration file see Section 6 3 4 6 2 updateSummary Var void updateSummaryVar const Vidx amp vIdx const NbrBox lt const UBA gentData gt 4 ubA gentDa taPtrNbrBox const Vector lt NbrBox lt REAL gt gt amp v gridPhiNbrBox const Vector lt NbrBox lt RE AL gt gt amp v_gridModelRealNbrBox cosnt Vector lt NbrBox lt 832 gt gt amp v_gridModelIntNbrBox V ector lt REAL gt amp v_realVar Vector lt S32 gt amp v_intVar Update the values of the summary variables The entire set of summary variables one variable per interface grid unit box for each attribute are reduced to get a summary value one summary value per attribute v dx indexes this unit box ubAgentDataPtrNbrBox holds pointers for UBAgent Data class instances for this unit box and its 26 neighboring unit boxes v_gridPhiNbrBox v gridModelRealNbrBox and v gridModelIntNbrBox pass the molecular concentrations the model specific REAL type attributes and the model specific S32 type attributes respectively for this unit box and its 26
42. detail sufficient to test their hypotheses translate the mathematical model into a computer program run the program on computers and visualize and analyze the simulation results While an ongoing challenge to this process is to gather and integrate sufficient biological knowl edge to build high fidelity models biologists have been studying the behavior of individual cells and subcellular organelles for a long time and have accumulated a great amount of knowledge about how cells and subcellular organelles behave The most natural way to build mathemati cal models based on this accumulated knowledge is to treat each cell or a subcellular organelle as a discrete simulation entity in combination with adoption of high resolution grids to model the environment However simulating these models often referred to as discrete agent based modeling is considered by some to be computationally intractable for the number of cells needed to study emergent behaviors of large living systems Kim et al 2007 Stamatakos 2010 More optimistic researchers have envisioned such simulations would be feasible on a supercomputer when a proper parallel implementation is used Jiao and Torquato 2011 We side with the latter researchers Biocellion is a software framework implemented by high performance computing HPC experts targeting large scale systems biology simulations and run ning on parallel computers ranging from laptops to supercomputers Efficient al
43. dx molecular species in the unit box gridRHS is initially set to the gridRHS value set inside the updatePDEBuf ferRHSLinear function for this molecular species Users may change the value based on model specifics 4 5 21 updatePDEBufferRHSTimeDependentSplitting void updatePDEBufferRHSTimeDependentSplitting const S32 pdeldx const VIdx amp startVIdx co nst Vidx amp pdeBufferBoxSize const Vector lt double gt amp v_gridPhi Vector lt double gt amp v_gridRHS Set the derivative values for the reaction term for the pdeldx PDE which is solved using the splitting scheme startVIdx points the low end corner of the PDE buffer unit box assuming the in terface grid spacing pdeBufferBoxSize is the size of a PDE buffer unit box assuming the interface grid spacing v_gridPhi holds molecular concentrations for the molecular species in this PDE the size of the vector is equal to the number of molecular species in this PDE and the first vector element holds the molecular concentration for the molecular species in this PDE with the smallest molecular species index This model routine updates v_gridRHS which holds the derivative values for the molecular species in this PDE 4 5 22 updatePDEBufferDirichletBC Val void updatePDEBufferDirichletBCVal const S32 elemidx const VReal amp startPos const VReal amp pdeBufferFaceSize const 32 dim const BOOL lowSide REAL amp bcVal Set the model specific D
44. e Vector lt REAL gt getModelRealArray void const Return a Vector array holding the entire set of model specific REAL type variable values void setModelReal const S32 idx const REAL val Set the value of the idx model specific REAL type variable to val void setModelRealArray const Vector lt REAL gt amp v val Set the entire set of model spe cific REAL type variables to v_val v_val s size should coincide with the number of model specific REAL type variables for this discrete agent type void incModelReal const 32 idx const REAL inc Increment the value of the idx model specific REAL type variable by inc e 532 getModellnt const S32 idx const Return the value of the idx model specific S32 type variable e Vector lt S32 gt getModellntArray void const Return a Vector array holding the entire set of model specific S32 type variable values void setModelInt const S32 idx const S32 val Set the value of the idx model specific S32 type variable to val e void setModelIntArray const Vector lt S32 gt amp v val Set the entire set of model specific S32 type variables to v val v val s size should coincide with the number of model specific 27 S32 type variables for this discrete agent type void incModelInt const S32 idx const S32 inc Increment the value of the idx model specific 32 type variable by inc 2 4 23 AgentMechIntrctData This class instance holds the sum of the forces on
45. e variables for this discrete agent type void incModelReal const S32 idx const REAL inc Increment the value of the idx model specific REAL type variable by inc e 32 getModellnt const S32 idx const Return the value of the idx model specific 32 type variable e Vector lt S32 gt getModelIntArray void const Return a Vector array holding the entire set of model specific S32 type variable values void setModelInt const S32 idx const 32 val Set the value of the idx model specific S32 type variable to val void setModelIntArray const Vector lt S32 gt amp v_val Set the entire set of model specific S32 type variables to v_val v_val s size should coincide with the number of model specific S32 type variables for this discrete agent type void incModelInt const S32 idx const S32 inc Increment the value of the idx model specific 32 type variable by inc 2 4 19 AgentJunctionInfo This class instance stores a list of junctions belong to this discrete agent and the state variables of the junction ends A unique ID is assigned to every discrete agent particle This ID is mainly 25 used internally but Biocellion allows users to read the ID if necessary If an agent is created in the previous baseline time step this class instance also holds the ID of the mother agent getPrevld returns the ID of the mother agent Otherwise getPrevld and getCurld return the same value Public non static memb
46. ecific S32 type attributes respectively about the extra cellular space for this box and the 26 neighboring boxes v_y stores the current values of the ODE system variables This model rou tine updates v_f which stores the derivatives for the variables in the ODE system The ODE sys tems for a discrete agent are updated before updateSpA gentState is called for the agent 4 3 3 updateSpAgentState void updateSpAgentState const Vidx amp vidx const AgentJunctionInfo amp junctionInfo const VReal amp vOffset const Vector lt NbrBox lt REAL gt gt amp v_gridPhiNbrBox const Vector lt NbrBox lt REAL gt gt amp v_gridModelRealNbrBox const Vector lt NbrBox lt S32 gt gt amp v_gridModelIntNbrBox SpAg entState amp state Update the state of this discrete agent v dx indexes the unit box this discrete agent is located in junctionInfo stores the list of junctions between this discrete agent and its neighboring dis crete agents vOffset is the offset from the center of the unit box this discrete agent belongs to v_gridPhiNbrBox v_gridModelRealNbrBox and v_gridModelIntNbrBox provide the information molecular concentrations model specific REAL type attributes and model specific S32 type at tributes respectively about the extra cellular space for this box and the 26 neighboring boxes This model routine updates state based on the provided information state initially holds the state of this discrete agent updateSpAgentState is calle
47. ed in vIdx1 indexes the unit box spAgentl is located in dir is a unit vector from the position of spAgent to the position of spAgent0 dist is the dis tance between the positions of the two discrete agents This model routine updates link endO and endl Set link to true to form a junction between the two discrete agents end is the junction end in the spAgent0 side and endl is the junction end in the spAgent side endO and end are irrelevant if link is set to false 4 4 2 computeForceSpAgent void compute ForceSpA gent const VIdx amp vIdx0 const SpAgent amp spAgent0 const Vector lt REAL gt amp v gridPhi0 const Vector lt REAL gt amp v_gridModelReal0 const Vector lt S32 gt amp v_gridModell nt0 const VIdx amp vIdx1 const SpAgent amp spAgentl const Vector lt REAL gt v_gridPhil const Ve ctor lt REAL gt v_gridModelReall const Vector lt S32 gt amp v_gridModellntl const VReal amp dir co nst REAL amp dist VReal amp force Update force between two discrete agents v dx0 indexes the unit box spAgent0 is located in v_gridPhi0 v_gridModelReal0 and v_gridModelIntO are molecular concentrations model specific REAL type variables and model specific S32 type variables for the unit box spAgent0 belongs 46 to respectively v dx indexes the unit box spAgent is located in v_gridPhil v_gridModel Reall and v gridModelInt are molecular concentrations model specific REAL type variables and model specific S32 type
48. ems in their entirety 4 1 2 Simulation Domain Biocellion supports a three dimensional rectangular simulation domain Biocellion assumes a rectangular grid the interface grid with a fixed grid spacing the interface grid spacing The interface grid spacing defines the resolution of extra cellular space Biocellion users set the interface grid spacing h Different computational modules communicate through this interface grid See Figures 1 4 and 1 5 For each grid aligned h x h x h box unit box Biocellion main tains a list of agents located within the box The interface grid spacing should be equal to or larger than the maximum direct physico mechanical interaction distance between any two discrete agents This assures that a discrete agent directly interacts only with discrete agents in the same box or the neighboring 26 boxes dX E H Bi y z VX a x y z X s x y z h interface grid spacing a x y 2 s x y z B x y z Figure 1 4 Biocellion asks users to set partial differential equation PDE parameters at the gran ularity of a unit box For reaction diffusion PDEs users set the decay rate A and the source term s for every unit box and the diffusion coefficient B for every face between two adjacent unit boxes sharing a face User provided model routines set PDE parameters for a unit box based on the state of the unit box and the states of the discrete agents located in the unit box For PDE parameters set fo
49. er functions e 64 getPrevid void const Return the ID of the mother agent if this discrete agent is newly created Otherwise this function and getCurld return the same value e 64 getCurld void const Return the ID of this agent e 32 getNumJunctions void const Return the number of junctions between this agent and other neighboring agents e S64 getOtherEndld const S32 idx const Return the ID of the discrete agent at the other end of the idx junction e const JunctionEnd amp getJunctionEndRef const S32 idx const Return a constant reference of the idx junction s end at this agent s side e JunctionEnd amp getJunctionEndRef const S32 idx Return a modifiable reference of the idx junction s end at this agent side e BOOL isLinked const AgentJunctionInfo amp otherInfo const Return true if there is a valid junction between this agent and another agent assuming that otherInfo is the AgentJunction Info type variable of the other discrete agent Return false if there is no junction between the two discrete agents e BOOL isLinked const AgentJunctionInfo amp otherInfo S32 amp idxThis S32 amp idxOther const Return true if there is a valid junction between this agent and another agent assuming that otherInfo is the AgentJunctionInfo type variable of the other discrete agent Return false if there is no junction between the two discrete agents This function returns the indices to access the junct
50. es cannot exceed the interface grid spacing v_spAgentState and v_spAgentDisp should have same size one element for one newly secreted discrete agent state initially holds the state of this discrete agent spAgentSecretionBySpAgent is called before updateSpA gent BirthDeath within a single baseline time step 44 4 3 5 updateSpAgentBirthDeath void updateSpAgentBirthDeath const Vidx amp vida const SpAgent amp spAgent const AgentMechInt rctData amp mechIntrctData const Vector lt NbrBox lt REAL gt gt amp v_gridPhiNbrBox const Vector lt NbrBox lt REAL gt gt amp v gridModelRealNbrBox const Vector lt NbrBox lt S32 gt gt amp v_gridModell ntNbrBox BOOL amp divide BOOL amp disappear Mark whether this discrete agent will divide or disappear v dx indexes the unit box this dis crete agent is located in spAgent provides the information about this discrete agent mechIntrct Data stores the temporary mechanical interaction state of this discrete agent v_gridPhiNbrBox v_gridModelRealNbrBox and v gridModelIntNbrBox provide the information molecular concen trations model specific REAL type attributes and model specific 32 type attributes respectively about the extra cellular space for this box and the 26 neighboring boxes This model routine updates divide and disappear Set divide to true if this discrete agent is going to divide di videSpA gent Section 4 3 7 will be called if divide is set to true Set disappear
51. et const S32 zOffset const Return true if the indexed box is in a valid interface region Return false if the box is located outside the simulation domain or in a PDE buffer region e BOOL getValidFlag const VIdx amp vOffset const Same to the above function Use vOffset instead of xOffset yOffset zOffset to index the unit box 28 const T amp getVal const S32 xOffset const S32 yOffset const S32 zOffset const Return a constant reference of the value associated with the indexed unit box This function should not be invoked for an invalid box const T amp getVal const Vidx amp vOffset const Same to the above function Use vOffset instead of xOffset yOffset zOffset to index the unit box void setVal const 32 xOffset const S32 yOffset const S32 zOffset const T val Set the variable associated with the indexed unit box to val This function should not be invoked for an invalid box void setVal const VIdx amp vOffset const T val Same to the above function Use vOffset instead of xOffset yOffset zOffset to index the unit box 2 4 26 IfGridBoxData This class instance holds a set of variables for a sub region of the simulation domain one value per one unit box in the partition This class instance is often used to pass a set of variables for a single partition or a single partition plus the valid ghost region of the partition T eent const VIdx amp vIdx const Return the value associated with the u
52. fault Every agent has an associated type Bio cellion users assign different numbers of model specific REAL type and S32 type variables or discrete agent attributes for different agent types e junctionEndType t 16 bit unsigned integer by default Every junction end Section 1 1 has an associated type Biocellion users assign different numbers of model specific REAL type and S32 type variables or junction end attributes for different junction end types e agentId_t 64 bit signed integer by default Every discrete agent has a unique ID e idx_t 16 bit signed integer by default Three idx t type variables are used to index a unit box the size of a unit box is h x h x h assuming that the grid spacing for the interface grid a The current version of Biocellion runs only on x86 compatible architectures This may change in the future 11 is h e ubAgentIdx t 16 bit unsigned integer by default This data type is used to access discrete agents in a single unit box We assume that the default sizes of the above data types are large enough for most practical appli cations If this does not hold users need to redefine the above data types and recompile the entire Biocellion framework and model libraries Advanced users can redefine the above data types to a smaller integer data type to save memory if the smaller data type is sufficient for their applications 2 3 Enumerations Biocellion has multiple enumerated data types which a
53. ghboring discrete agents Figure 1 2 Each discrete agent has a set of junction ends one for each junction and each junction end has associated state variables Biocellion asks users to set the number of junction end types and the number of attributes or state variables for each junction end type Every junction automatically breaks if the distance between a pair of discrete agents exceeds the user provided maximum direct interaction distance for the pair Figure 1 2 A junction between two discrete agents Each discrete agent has a set of variables storing its temporary mechanical interaction state one 3 default variable is force users can add extra variables based on model requirements Temporary mechanical interaction state variables are valid within a single baseline time step time steps are explained in Section 1 3 Temporary mechanical interaction state variables are updated while evaluating physico mechanical interactions between discrete agent pairs These variables are used to set discrete agent displacement and to update agent state based on model specific rules The surrounding environment influences biological system behavior The surrounding environ ment affects the transport of molecules in the extra cellular space and promotes and limits cell movement e g chemotaxis cells moving towards either the upward of downward gradient of molecular concentrations promotes cell movement and soil aggregate structure constrains
54. gorithms careful performance tuning and exploitation of parallelism at multiple levels increases computing capabil ity well beyond that of general simulation frameworks enabling simulation of high fidelity models at large scale often several orders of magnitude faster than on other simulation platforms Biocel lion relieves computational biologists of the burdens posed by performance and scalability freeing them to focus instead on the modeling task itself Biological systems are diverse and highly complex and there are not yet any standards for math ematical modeling flexible enough to express the expansive biological knowledge integrated into component models Mathematical models depend on the specifics of the biological systems ques tions being answered required level of accuracy availability of biological knowledge and model developer s intellectual capability and experience Software designed for a single mathematical model will have only a short life span Each model is unique in its design While the need for generality in modeling is therefore paramount we have noticed that many multicellular system models exhibit similar computational challenges We have studied a wide range of biological systems and analyzed numerous biological system models to separate individual model specifics from common computational challenges We have abstracted out and integrated into the Biocel lion framework what is common to a wide spectrum of life sy
55. h the same variable name If users need to set shared global data across the entire set of MPI processes update the data using the updateGlob alData model routine Section 4 2 14 which is invoked exactly once by just one of the MPI processes the MPI process having rank 0 The data is copied to the remaining MPI processes by the Biocellion core framework The data cannot be modified outside the updateGlobalData function and users can access the data through getGlobalDataRef Section 3 1 4 Section 4 1 summarizes model routine invocation timings and the remaining of this Section ex plains the model routines 4 1 Model Routine Invocation Timings Figures 4 1 4 2 4 3 and 4 4 summarizes model routine invocation timings for the initialization main computation PDE computation and termination respectively Symbols after model rou tine names indicate the model routine invocation granularity S indicates that the model routine is called only once per simulation even when there are more than one MPI processes M indicates that the model routine is called once per MPI process If there is only one multi threaded pro cess S and M are identical P indicates that the model routine is called once per partition IB indicates that the model routine is called once per interface grid unit box IF indicates that the model routine is called once per relevant unit box face e g updatelfGridBetaDomainBdry is called only for the unit box faces at t
56. he domain boundary PB indicates that the model routine is called once per PDE buffer grid unit box and PF indicates that the model routine is called once per relevant PDE buffer grid unit box face A indicates that the model routine is called once per discrete agent and AP indicates that the model routine is called once per discrete agent pair within the maximum direct interaction distance 37 updatelfGridSpacing S UpdateOptModelRoutineCallinfo S updatgeDomainBdryType S updatePDEBufferBdryType S updateTimeStepInfo S updateSyncMethod S updateSpAgentinfo S updateJunctionEndInfo S 9 updatePDElnfo S 10 updatelfGridModelVarInfo S 11 updateRNGinfo S 12 updateFileOutputinfo S 13 updateSummaryOutputinfo S Start initialization Parse the simulation configuration file SoS Provide model information Setup global initialization data Perform model specific initialization ON AMAWNE 1 updateGlobalData S 1 init M Initialize using check point data Set PDE buffer regions Set habitable regions initialize agent and grid data End initialization Figure 4 1 Initialization No Read check point data 1 setPDEBuffer P setHabitable IB addSpAgents P initJunctionSpAgent AP initlfGridVar IB initPDEBufferPhi PB Initialize using check point data um P DH 4 2 Mode
57. he expo nential distribution is selected and P if the gamma distribution is selected e REAL param2 Relevant only when the gamma distribution is selected The displacement value for the gamma distribution 2 4 15 FileOutputInfo This class instance holds the information related to file output of simulation results Public non static member variables e BOOL particleOutput If set to true Biocellion generates output files for discrete agent data 22 e S32 particleNumExtraOutputVars This variable sets the number of extra output variables in addition to the default radius and color variables for particles representing discrete agents e Vector lt BOOL gt v_gridPhiOutput Set to true to generate output files for each molecular species in the PDE e Vector lt BOOL gt v_gridPhiOutputDivideByKappa Set to true to output molecular concen tration values divided by k values Section 4 5 3 for each molecular species in the PDE 2 4 16 SummaryOutputInfo This class instance holds the information related to summary output of simulation results Public non static member variables e Vector lt string gt v realName Each vector element sets the name of the corresponding sum mary attribute e Vector lt summary type e gt v realType Each vector element sets the reduction method to print summary for the corresponding summary attribute 2 4 17 SpAgentState This class instance holds variables describing the state of a discrete agent and p
58. ibmodel model yeast pattering model define h should be set to 1 before compiling the model code To adopt the 20 um grid spacing H GR D_SPAC NG DM CRO_METER needs to be set to 0 57 6 0 Biocellion Installation and Execution 6 1 System Requirements Biocellion runs on a spectrum of parallel computers ranging from multicore PCs and workstations to cluster computers Cloud computers and supercomputers users do not need to change their model code to run on different platforms The current version of Biocellion runs only on x86 compatible systems a great majority of computers with Intel or AMD CPUs are x86 compatible Biocellion runs on 64 bit Linux operating systems Compiling Biocellion model code requires a C compiler we have tested using the GNU gcc compiler and the Intel icc compiler Biocellion uses the Intel Thread Building Blocks library for multi threading and the library is freely available from the thread building blocks homepage http threadingbuildingblocks org Biocellion is pre installed in Amazon EC2 Section 6 6 6 2 Compiling Biocellion 6 2 1 Compiling Model Code Biocellion users first need to update directory paths in SBTOCELLION_ROOT Makefile com mon BIOCELLION_ROOT should point the Biocellion root directory REAL in SBIOCELLION ROOT Makefile common is set to double by default Users can set REAL to float to use single precision f
59. ies in this PDE the size of the vector is equal to the number of molecular species in this PDE and the first vector element holds the molecular concentration for the molecular species in this PDE with the smallest molecular species index v_gridModelReal and v gridModelInt are the model specific REAL type variables and model specific S32 type variables for this unit box respectively This model routine updates v_gridRHS which holds the derivative values for the molecular species in this PDE 4 5 11 updatelfGridAMRTags void updatelfGridAMRTags const Vidx amp vIdx const NbrBox lt const UBAgentData gt amp ubA gent DataPtrNbrBox const Vector lt NbrBox lt REAL gt gt amp v_gridPhiNbrBox const Vector lt NbrBox lt REAL gt gt amp v_gridModelRealNbrBox const Vector lt NbrBox lt S32 gt gt amp v_gridModelIntNbrBox Vector lt S32 gt amp v finestLevel Set the desired adaptive mesh refinement AMR level for this unit box v dx indexes this unit box ubAgentDataPtrNbrBox holds pointers for UBAgentData class instances for this unit box and its 26 neighboring unit boxes v gridPhiNbrBox v_gridModelRealNbrBox and v_gridMod elIntNbrBox are the molecular concentrations model specific REAL type attributes and model specific 32 type attributes for this unit box and its 26 neighboring unit boxes respectively This model routine updates v_finestLevel which holds the desired finest AMR level for this unit box for the PDEs in this simulati
60. inVerbosity msg Print msg to the standard output if the standard output verbosity set inside the simulation configu ration file see Section 6 3 is equal to or larger than minVerbosity 34 3 4 3 ERROR ERROR msg Print msg to the standard output and abort the program 35 4 0 Biocellion Model Routines Biocellion asks users to provide their model specifics by filling the function body of a set of prede fined model routines Users write only sequential code but the Biocellion framework can invoke a single model routine multiple times in parallel Users are disallowed to use any functions that are not thread safe e g the C rand function is not thread safe use getModelRand Section 3 2 1 instead We also strongly discourage users to access global variables inside a model routine with few ex ceptions Reading constant global variables is acceptable If global variables are initialized at the beginning of the simulation in the init model routine Section 4 2 15 reading the initialized vari ables is also acceptable with one caveat If multiple MPI processes are created to run Biocellion on a system with multiple shared memory nodes e g cluster computers and supercomputers the init model routine is called once in every MPI process Model routines executed in different MPI processes access different instances of a global variable initialized in different init invoca tions even when the routines access a variable wit
61. ion ends in addition idxThis holds the index for the junction end at this discrete agent s side and idxOther holds the index for the junction end at the other discrete agent s side 2 4 20 SpAgent This class instance holds the position offset within the interface grid unit box this agent belongs to and the state of a cell or a non cellular biological entity mapped to a single particle Public non static member variables e SpAgentState state A variable holding the state of this agent e AgentJunctionInfo junctionInfo A variable holding information about the entire set of 26 junctions between this agent and neighboring agents e VReal vOffset This variable holds the offset of this discrete agent from the center of the interface grid unit box this discrete agent belongs to 2 4 21 UBAgentData This class instance holds a list of discrete agents in an interface grid unit box e Vector lt SpAgent gt v_spAgent A collection of discrete agents in this unit box 2 4 22 ExtraMechIntrctData This class instance holds the extra temporary mechanical interaction state variables of a discrete agent The entire set of extra temporary mechanical interaction state variables are stored in an internal linear array Biocellion does not allow direct access to the array which can be error prone Public non static member functions e REAL getModelReal const S32 idx const Return the value of the idi model specific REAL type variable
62. ionSize e Vidz getldx1DTo3D const S64 idx const idx_t size Perform the reverse operation of getldx3DTo1D const Vidx amp vidx const idx_t size 2 4 3 VReal VReal is a container for three REAL type variables Public static member variables e VReal ZERO A constant three dimensional zero vector i e 0 0 0 0 0 0 e VReal UNIT A constant three dimensional unit vector i e 1 0 1 0 1 0 e VRealA BASIS 3 A set of three dimensional basis vectors ABASIS 0 1 0 0 0 0 0 A BASIS 1 0 0 1 0 0 0 and A BASIS 2 0 0 0 0 1 0 Public non static member functions e x Arithmetic operators e idx An access operator returns the idx element of the internal array holding three REAL type variables idx should be 0 1 or 2 16 2 4 4 ODENetInfo This class instance holds the information required to set up a solver for a deterministic ordinary different equation ODE system to model a discrete agent regulatory network Public non static member variables e 32 numVars Set the number of variables in the ODE system e ode_stiff_e stiff Set the stiffness of the ODE system Biocellion uses different numerical algorithms based on this information e REAL h Set the initial step size to h x state and grid time step size 0 0 lt h lt 1 0 e REAL hm Set the minimum time step size to hm x state and grid time step size 0 0 lt hm lt h e REAL epsilon Set the relative error tolerance
63. irichlet boundary value for the elemIdx molecular species for a PDE buffer unit box face at the domain boundary This model routine is valid only when BC TYPE DIRICHLET MODEL is selected startPos points the low end corner of the unit box face pde BufferFaceSize is the size of a PDE buffer unit box face The unit box face is orthogonal to a unit vector in the x if dim is 0 y if dim is 1 or z if dim is 2 direction The face is at the low side if lowSide is true or the high side if lowSide is false of the simulation domain This model routine updates bcVal 4 5 23 updatePDEBufferNeumannBCVal void updatePDEBufferNeumannBCVal const 32 elemldx const VReal amp startPos const VReal amp pdeBufferFaceSize const S32 dim const BOOL lowSide REAL amp bcVal Set model specific Neumann boundary values for the elemidx molecular species for a PDE buffer unit box face at the domain boundary This model routine is valid only when BC TYPE NEUMANN MODEL is selected startPos points the low end corenr of the unit box face pdeBuffer FaceSize is the size of a PDE buffer unit box face The unit box face is orthogonal to a unit vector 54 in the x if dim is 0 y Gf dim is 1 or z if dim is 2 direction The face is at the low side if lowSide is true or the high side if lowSide is false of the simulation domain This model routine updates bcVal 4 6 Simulation Output 4 6 1 updateSpAgentOutput void up
64. l Configuration 4 2 1 updatelfGridSpacing void updatelfGridSpacing REAL amp ifGridSpacing Set the interface grid spacing This model routine sets ifGridSpacing to the desired interface grid spacing 4 2 2 updateOptModelRoutineCallInfo void updateOptModelRoutineCallInfo OptModelRoutineCallInfo amp callInfo Control the invocation of optional model routines This model routine updates callInfo See Sec tion 2 4 6 about the OptModelRoutineCallInfo data type 4 2 3 updateDomainBdryType void updateDomainBdryType domain_bdry_type_e a domainBdryType 3 Set the domain boundary types at the low and high sides in the x y and z directions This model routine updates a_domainBdryType 3 See the explanation about domain_bdry_type_e in Sec tion 2 3 for the supported domain boundary types 38 1 updateSpAgentOutput A Start main Startmainloop 2 updateSummaryVar IB updatelfGridKappa IB updatelfGridAlpha IB updatelfGridBetalnifRegion IF updatelfGridBetaGridBufferBdry IF updatelfGridBetaDomainBdry IF NONN adjustifGridRHSTimeDependentLinear 18 8 updatelfGridRHSTimeDependentSplit i ting IB Ter Si 9 updatelfGridDirichletBCVal IF baseline time 10 updatelfGridNeumannBCVal IF steps 11 updatePDEBufferKappa PB 12 updatePDEBufferAlpha PB Yes 13 updatePDEBufferBetalnGridBufferRe gion PF Inject cells 14 updatePDEBufferBetaDomainBdry PF iterl lt 15 updatePDEBuffer
65. l equation ODE and partial differential equation PDE solver outputs to assure that output values are equal to or larger than the user specified minimum valid value Sections 2 4 4 and 2 4 9 If PDE solver output values contain a value smaller than the specified minimum value Biocellion prints the PDE parameters and the molecular concentrations of the unit box with the erroneous value and the neighboring boxes in the x y and z directions at the beginning of the PDE time step along with the erroneous value after the PDE time step which produced the erroneous value Biocellion aborts after printing the error message 6 8 Visualization Using Paraview We briefly summarize visualization instructions assuming Paraview 3 14 1 See the Paraview homepage and manual for advanced visualization 6 8 1 Visualizing Discrete Agents 1 Click the File menu and select the Open item 2 Open a set of files in the parallel VTK unstructured grid format pvtu 63 10 6 8 2 Click the Properties tab at the bottom left side and click the Apply icon Select the Filters menu the Alphabetical sub menu and the Glyph item Select the Properties tab Set Scalars to radius Glyph Type to Sphere Radius to 1 0 Scale Mode to Scalar Set Scale Factor to 1 0 and Maximum Number of Points to a value larger than the number of discrete agents Click Apply
66. l specific S32 type variables for this discrete agent type void incModelInt const S32 idx const S32 inc Increment the value of the idx model specific 32 type variable by inc 2 4 18 JunctionEnd This class instance holds variables describing the state of a junction end and provides interface functions to access the variables The entire set of junction end state variables are stored in an 24 internal linear array Biocellion does not allow direct access to the array which can be error prone Public non static member functions e junctionEndType_t getType void const Return the type of the junction end void setType const junctionEndType_t type Set the type of the junction end to type This function also resizes and resets the internal linear array Users need to properly initialize junction end state variables after this function is invoked e REAL getModelReal const S32 idx const Return the value of the idx model specific REAL type variable e Vector lt REAL gt getModelRealArray void const Return a Vector array holding the entire set of model specific REAL type variable values void setModelReal const S32 idx const REAL val Set the value of the idx model specific REAL type variable to val void setModelRealArray const Vector lt REAL gt amp v val Set the entire set of model spe cific REAL type variables to v val v val s size should coincide with the number of model specific REAL typ
67. ler than the interface grid spacing The maximum direct interaction distance between two discrete agents having different types is computed by averaging the dMax values for the two types In computing direct interactions between discrete agent pairs every pair within the maximum interaction distance is enumerated e S32 numStateModelReals Set the number of model specific REAL type state variables for this discrete agent type e 32 numsStateModellnts Set the number of model specific S32 type state variables for this discrete agent type e S32 numExtraMechIntrctModelReals Set the number of model specific REAL type tem porary mechanical interaction state variables for this discrete agent type e S32 numExtraMechIntrctModelInts Set the number of model specific S32 type temporary 18 mechanical interaction state variables for this discrete agent type 2 4 8 JunctionEndInfo This class instance holds the model specifics about a junction end type Public non static member variables e S32 numModelReals Set the number of model specific REAL type variables for this junc tion end type e S32 numModellnts Set the number of model specific 32 type variables for this junction end type 2 4 9 GridPhilnfo This class instance holds the model specifics about the concentration values of a molecular species which are mainly updated by solving PDEs Public non static member variables e S32 elemldx Set the unique index for this molecular s
68. lfRegion IF updatelfGridBetaGridBufferBdry IF updatelfGridBetaDomainBdry IF updatelfGridRHSLinear IB updatePDEBufferKappa PB updatePDEBufferAlpha PB updatePDEBufferBetalnGridBufferRe gion PF 10 updatePDEBufferBetaDomainBdry PF pdatePDEBufferRHSLinear PB Start solving a PDE Set PDE parameters 070 ACUTE Nees 1 adjustlfGridRHSTimeDependentLinear IB 2 adjustPDEBufferRHSTimeDependentLinear PB type steady state 1 updatelfGridRHSTimeDependentSplitting IB reaction 2 updatePDEBufferRHSTimeDependentSplitting PB diffusion linear Update boundary condition Solve an elliptic linear PDE type time dependent reaction diffusion linear No type time dependent reaction diffusion splitting iterO lt PDE time steps Solve a system of ODEs iterO lt PDE time steps Update boundary condition Adjust source parameters Ll A single step a parabolic PDE splitting PDE diffusion time steps Update boundary condition Single step a parabolic PDE updatelfGridDirichletBCVal IF updatelfGridNeumannBCVal IF updatePDEBufferDirichletBCVal PF updatePDEBufferNeumannBCvVal PF End solving a PDE PDA Figure 4 3 Solving a PDE Start termination Perform model specific termination End initialization Figure 4 4
69. lnlfRegion void updatelfGridBetalnlfRegion const 32 elemldx const S32 dim const VIdx amp vIdx0 const VI dx amp vidxl const UBAgentData amp ubAgentData0 const UBAgentData amp ubAgentDatal const Ve ctor lt REAL gt v_gridPhi0 const Vector lt REAL gt v gridPhil const Vector lt REAL gt v grid ModelReal0 const Vector lt REAL gt amp v_gridModelReall const Vector lt S32 gt amp v gridModelIntO const Vector lt S32 gt amp v_gridModellntl REAL amp gridBeta Set the diffusion coefficient for the elemidx molecular species for a unit box face inside the interface region The unit box face is orthogonal to a unit vector in the x if dim is 0 y if dim is 1 or z Gif dim is 2 direction vidx0 indexes the unit box at the low side of the face and vidx1 indexes the unit box at the high side of the face ubAgentData0 and ubAgentDatal hold discrete agents at the low and high side unit boxes respectively v gridPhi0 and v_gridPhil are molecular concentrations for the low and high side unit boxes respectively v_gridModelReal0 and v_grid ModelReall are the model specific REAL type variables for the unit boxes at the low and high sides v_gridModellnt0 and v_gridModelInt1 are the model specific 32 type variables for the unit boxes at the low and high sides respectively This model routine updates gridBeta 4 5 6 updatelfGridBetaPDEBufferBdry void updatelfGridBetaPDEBufferBdry const 32 elemIdx const S32 dim const VIdx amp vida
70. loating point arithmetic instead of double precision floating point arithmetic MPAGENT should be set to no and SPAGENT should be set to yes in the current version Users can set DEBUG to yes for debugging support to use GDB or no for high performance Users may update several variables in SBIOCELLION ROOT Makefile model to use differ ent compilers or different compilation flags the default is to use GNU gcc We strongly encour age users to enable checks on model routine outputs and input arguments of Biocellion support functions by setting CHECK FLAG to DENABLE CHECK 1 at the early stage of model devel opment Once the model code is well verified users may disable the checks set CHECK FLAG to DENABLE CHECK 0 to expedite simulation Biocellion builds a model library using the model code under SBIOCELLION_ROOT 1libmo del model To test using one of the two sample models users create a soft link e g under SBIOCELLION ROOT libmodel type In s model yeast patterning updated model and type make under SBIOCELLION ROOT 1libmodel Once compilation completes 1ibmodel so will be created under SBIOCELLION_ROOT 1libmodel interface 6 2 2 Compiling the Biocellion Core Framework Advanced users may want to recompile the Biocellion core framework to use a different compiler with different optimization flags or to link the Biocellion
71. lobalDataRef const Vector lt U8 gt getGlobalDataRef void Return a constant reference to the user global data Biocellion invokes the model routine update GlobalData see Section 4 2 14 to set the global data this data can be accessed in any model routine once the data is set during the simulation initialization process This function returns a constant reference to the data 3 1 5 getRecentSummaryRealVal REAL getRecentSummaryRealVal const S32 elemidx Return the most recent value of the elemIdx REAL type summary variable 31 3 1 6 getRecentSummaryIntVal REAL getRecentSummaryIntVal const S32 elemidx Return the most recent value of the elemIdx 32 type summary variable 3 1 7 getDomainSize idx 1 getDomainSize const S32 dim Return the domain size in the x if dim is 0 y Gif dim is 1 or z if dim is 2 direction The domain size is set inside the simulation configuration file see Section 6 3 3 1 8 getPartSize idx_t getPartSize void Return the partition size set inside the simulation configuration file see Section 6 3 3 1 9 getSimInitType sim_init_type_e getSimInitType void Return the simulation initialization type set inside the simulation configuration file see Sec tion 6 3 See the explanation about sim init type e in Section 2 3 for the valid simulation ini tialization types 3 1 10 getParticleNumExtraOutputVars S32 getParticleNumExtraOutputVars void Return the number of ex
72. m const VIdx amp s tartVIdx0 const Vidx amp startVIdx1 const VIdx amp pdeBufferBoxSize REAL amp gridBeta Set the diffusion coefficient for the elemldx molecular species for a PDE buffer grid unit box face inside a PDE buffer region The unit box face is orthogonal to a unit vector in the x if dim is 0 y Gif dim is 1 or z if dim is 2 direction startVIdx0 points the low end of the unit box at the low side of the face and startVldx1 points the low end of the unit box at the high side of the face assuming the interface grid spacing pdeBufferBoxSize is the size of a PDE buffer unit box assuming the interface grid spacing This model routine updates gridBeta 4 5 18 updatePDEBufferBetaDomainBdry void updatePDEBufferBetaDomainBdry const S32 elemIdx const 32 dim const VIdx amp startVId x const Vidx amp pdeBufferBoxSize REAL amp gridBeta Set the diffusion coefficient for the elemIdx molecular species for a PDE buffer unit box face at the domain boundary The unit box face is orthogonal to a unit vector in the x if dim is 0 y if dim is 1 or z if dim is 2 direction startVIdx points the low end corner of the PDE buffer unit box in the PDE buffer partition assuming the interface grid spacing pdeBufferBoxSize is the size of a PDE buffer unit box assuming the interface grid spacing This model routine updates gridBeta 4 5 19 updatePDEBufferRHSLinear void updatePDEBufferRHSLinear const 32 elemldx const Vidx amp sta
73. mary output Added the updateSummaryOutputInfo model routine Section 4 2 13 and the Summa ryOutputInfo class Section 2 4 15 Removed the printSummary and summaryType member variables from the GridPhilnfo and IfGridModel VarInfo classes Sections 2 4 9 and 2 4 13 respectively updateSum maryOutputinfo controls summary output Added the updateSummaryVar model routine Section 4 6 2 e Modify the Biocellion API related to setting options for partial differential equation PDE multigrid and ordinary differential equation ODE solvers Added the MGSolveInfo class Section 2 4 10 Added the mgSolveInfo member variable to the PDEInfo class Section 2 4 12 This variable controls options for the PDE multigrid solver Added the epsilon and threshold member variables to the ODENetInfo class Sec tion 2 4 4 Added the odeEpsilon and odeThreshold member variables to the SplittingInfo class Section 2 4 11 Removed the mg_elliptic mg_parablic and ode_solve elements Section 6 3 2 e Modified the syntax of the computeForceSpAgent Section 4 4 2 and computeExtraMechIn trctSpAgent Section 4 4 3 model routines The grid attributes of the unit boxes containing an interacting discrete agent pair are passed in addition Modified the syntax of the updatePDEBufferDirichletBC Val Section 4 5 22 and updatePDE BufferNeumannBCVal Section 4 5 23 model routines a_gridPhi is removed Added the update FileOutputInfo updateSumm
74. n about a single model specific attribute See Section 2 4 13 for additional details about the IfGridModel VarInfo class 4 2 11 updateRNGInfo void updateRNGInfo Vector lt RNGInfo gt amp v_rngInfo Provide the information about the random number generators used in the model This model routine updates v_rnglnfo The size of v_rngInfo coincides with the number of random number generators Each vector element provides the information about a single random number genera tor See Section 2 4 14 for additional details about the RNGInfo class 4 2 12 updateFileOutputInfo void update File OutputInfo FileOutputInfo amp fileOutputInfo Provide the information about file output of simulation results This model routine updates file OutputInfo See Section 2 4 15 for additional details about the FileOutputInfo class 4 2 13 updateSummaryOutputInfo void updateSummaryOutputinfo SummaryOutputInfo amp summaryOutputinfo Provide the information about summary output of simulation results This model routine updates summaryOutputinfo See Section 2 4 16 for additional details about the Summary OutputInfo class 4 2 14 updateGlobalData void updateGlobalData Vector lt U8 gt amp v_globalData Update the global data shared by every MPI process Once v_globalData is updated in this model routine the data is copied to every MPI process v_globalData cannot be updated elsewhere getGlobalDataRef Section 3 1 4 returns a constant reference to the
75. neighboring unit boxes This model routine updates summary values for REAL type attributes v_realVal and S32 type attributes v_intVal 55 5 0 Sample Implementations We provide two sample implementations porting Momeni et al s yeast patterning model Mo meni et al 2013 and Resat et al s bacteria on soil aggregate model Resat et al 2011 See the cited papers for additional model details Kang etal Kang et al TBP also portray the porting process of the yeast patterning model and we refer interested readers to the article The yeast patterning model code is located under B yeast pattering and the bacteria on soil aggregate model is located under SBIOCELL ON_ROOT 1libmodel model soil bacteria See SB k main yeast patterning 5um xml and B OCELL ON_ROOT libmodel model OCELL ON_ROOT framewor OCELL soil aggregate model is in SBIOCELL ml M ON_ROOT framework main yeast patterning 20um xml for simulation configuration files for the yeast patterning model with 5 um and 20 um interface grid spacings respectively see Section 6 3 for details about simulation configuration files A sample simulation configuration file for the bacteria on ON ROOT framework main soil bacteria x CRO METER in SB OCELL To adopt the 5 um grid spacing for the yeast patterning model H GRID SPACING 5 ON_ROOT 1
76. ng REAL epsilon The PDE multigrid solver assumes solution convergence if the residual norm is equal to or smaller than the original norm multiplied by epsilon e REAL hang Set the minimum required change in two consecutive V cycles or W cycles If the decrease in the norm relative to the norm in the previous iteration is less than hang the PDE multigrid solver exits prematurely and prints a warning REAL normThreshold The PDE multigrid solver assumes solution convergence if the residual norm is equal to or smaller than normThreshold 2 4 11 SplittingInfo This class instance holds the information required to solve a PDE using the splitting scheme A PDE can have multiple molecular species if the splitting scheme is used Public non static member variables 20 e Vector lt S32 gt v_diffusionTimeSteps Set the number of diffusion time step per PDE time step for each molecular species in the PDE e ode stiff e odeStiff Set the stiffness of the ODE system for the PDE reaction part Biocel lion uses different numerical algorithms based on this information e REAL odeH Set the initial ODE time step size to odeH x PDE time step size 0 0 lt odeH lt 1 0 e REAL odeHm Set the minimum ODE time step size to odeHm x PDE time step size 0 0 lt odeHm lt odeH e REAL odeEpsilon Set the relative error tolerance The ODE solver cannot ensure re quired accuracy if epsilon is smaller than le 9 e REAL odeThreshold If
77. ng partial differential equation PDE parameters Aunmtbozandits26neighboringbozes Simulation domain Tana Ka WAA ki WALIA Mapping a cell to multiple discrete agents Initialization o nie bot RIS GIS bw NTE hk OR E ed Main loop Solving a PDE ico Gs Se SS be SE ae ee Oe Se SS Termination xiv 38 39 40 40 1 0 Introduction Many biological phenomena arise as the product of complex interactions in living systems of cells These emergent behaviors are evident only once these systems grow to millions to trillions of cells Understanding and predicting emergent behaviors is crucial to solving important problems in the areas of energy Ha et al 2011 Majors et al 2008 Melnicki et al 2013 environment Daly 2000 Singh et al 2011 and healthcare Hall Stoodley et al 2004 Weinberg 2007 Rubinstein et al 2013 Computational biologists applying the scientific method iterate using modeling and simu lation refining their models until simulation results match those of experiments Biocellion pro vides the fast simulation of living systems needed to facilitate research into emergent behaviors Computer simulation of biological systems requires multiple steps Computational biologists must acquire biological knowledge or generate hypotheses on various aspects of biological system be haviors build a mathematical model based on this biological knowledge with
78. nit box at v dx The the sub region for this class instance should include v dx void set const Vidx amp vIdx const T val Set the value associated with the unit box at vidx to val The the sub region for this class instance should include v dx VIdx smallEnd void const Return the grid index for the unit box at the low end corner of the sub region for this class instance VIdx size void const Return the size of the sub region for this class instance idx_t size const S32 dim const Return the size of the sub region for this class instance in either x if dim is 0 y if dim is 1 or z if dim is 2 direction a The valid ghost region for a partition box can be obtained if we increase the size of the box by 1 in the x y and z directions the box size in each direction is increased by 2 and exclude the unit boxes located in the original box in the PDE buffer region and outside the simulation domain 29 3 0 Biocellion Support Functions 3 1 Simulation Instance Information 3 1 1 getCurBaselineTimeStep S32 getCurBaselineTimeStep void Return the current baseline time step 3 1 2 getCurStateAndGridTimeStep S32 getCurStateAndGridTimeStep void Return the current state and grid time step within the current baseline time step 3 1 3 getCurPDETimeStep S32 getCurPDETimeStep const S32 pdeldx Return the current PDE time step of the pdeldx PDE within the current state and grid time step 3 1 4 getG
79. of a simulated biological system Discrete agents may interact with each other directly and indirectly Direct interaction includes short range physico mechanical interaction between neighboring discrete agents such as cell cell adhesion and shoving or immune cells engulfing other discrete agents phagocytosis We assume the radius of direct interaction is relatively short e g within a small multiple of cell diameter d d d radius r height h direction d d d radius r Figure 1 1 Discrete agent shapes Indirect interactions result from the secretion of diffusible molecules and cells uptaking the se creted diffusible molecules Cell surface molecules receptors also bind to signaling molecules secreted by other cells Biocellion users set the maximum direct physico mechanical interaction distance for each discrete agent type The maximum interaction distance between two discrete agents with different types is computed by averaging the two maximum direct interaction distances for the two types Biocel lion invokes model routines that evaluate physico mechanical interaction between a cell pair only for discrete agent pairs within the maximum direct interaction distance Physico mechanical interactions between cell pairs often depend on past interactions Once two cells form an adherens junction pulling two cells apart requires stronger force Biocellion users can explicitly form and break a junction between a pair of nei
80. on index the vector using a PDE index not a molecular species index 4 5 12 updatelfGridDirichletBCVal void updatelfGridDirchletBCVal const 32 elemldx const VReal amp pos const S32 dim const BO OL lowSide const REAL a_gridPhi 3 const Vector lt Vector lt REAL gt gt amp vv_gridModelReal co nst Vector lt Vector lt S32 gt gt amp vv_gridModellnt REAL amp bcVal Set the model specific Dirichlet boundary value for the elemIdx molecular species for a unit box face at the domain boundary This model routine is valid only when BC TYPE DIRICHLET MODEL is selected pos locates the center of the unit box face The unit box face is orthogonal to a unit vector in the x if dim is 0 y if dim is 1 or z if dim is 2 direction The face is at the low side if lowSide is true or the high side if lowSide is false of the simulation domain a_gridPhi holds molecular concentrations for the elemIdx molecular species for the three interface grid unit boxes closest to the boundary unit box face along the direction orthogonal to the face vv_grid ModelReal and vv gridModelInt hold model specific REAL and S32 type variables for the three boxes closest to the boundary unit box face e g vv_gridModelReal elemldx 0 holds the value for the elemIdx model specific REAL type variable for the unit box closest to the boundary face 51 This model routine updates bcVal 4 5 13 updatelfGridNeumannBC Val void updatelfGridNeum
81. on the entire simulation domain pdeBufferBoxSize is the size of a PDE buffer unit box assuming the interface grid spacing This model routine updates v_gridPhi 4 5 15 updatePDEBufferKappa void updatePDEBufferKappa const S32 pdeldx const Vidx amp startVIdx const VIdx amp pdeBufferB oxSize REAL amp gridKappa Set the ratio of the available volume used in computing diffusion flux for a unit PDE buffer box of the pdeldx PDE startVIdx points the low end corner of the PDE buffer unit box assum ing the interface grid spacing pdeBufferBoxSize is the size of a PDE buffer unit box assuming the interface grid spacing This model routine updates gridKappa 0 0 lt gridKappa lt 1 0 Set gridKappa to 1 0 to assume that 100 of the volume is free gridKappa is assumed to be 1 0 outside the simulation domain 4 5 16 updatePDEBufferAlpha void updatePDEBufferAlpha const S32 elemldx const Vidx amp startVIdx const Vidx amp pdeBufferB oxSize REAL amp gridAlpha Set the decay rate for the elemidx molecular species startVIdx points the low end corner of 52 the PDE buffer unit box assuming the interface grid spacing pdeBufferBoxSize is the size of a PDE buffer unit box assuming the interface grid spacing This model routine updates gridAlpha 1 0 lt gridAlpha lt 0 0 Set gridAlpha to 0 0 to ignore decay 4 5 17 updatePDEBufferBetaInPDEBufferRegion void updatePDEBufferBetaInPDEBufferRegion const 32 elemidx const S32 di
82. ormation about a molecular species in this PDE If the splitting scheme is used there can be more than one molecular species in a single PDE 2 4 13 IfGridModelVarInfo Biocellion allows users to assign a fixed number of model specific attributes to the interface grid or a fixed number of REAL and S32 type state variables to every unit box in the interface grid These variables are updated by model specific rules in the model routine called at the beginning and the end of every state and grid time step This class instance provides the information about a model specific attribute Public non static member variables e std string name Set the name of the attribute 2 4 14 RNGInfo Biocellion provides random number generators with the uniform Gaussian exponential and gamma distributions We assume that the probability density functions for the exponential and gamma EE 1 distributions are B x e P and B x T o x x x e BX respectively Public non static member variables e rmg_type_e type This variable sets the random number generator type e REAL param0 The minimum value if the uniform distribution is selected the average value if the Gaussian distribution is selected B if the exponential distribution is selected and a if the gamma distribution is selected e REAL param The maximum value if the uniform distribution is selected the standard deviation value if the Gaussian distribution is selected the displacement value if t
83. ox radius should not exceed one half of the interface grid spacing ubVIdxOffset idx idx 0 1 or 2 should be 1 0 or 1 3 3 2 changePosFormat2LvTolLv void changePosFormat2LvTolLv const VIdx amp vIdx const VReal amp vOffset VReal amp pos Particle position can be represented in a single level with a single VReal varaible storing the offset from the low end corner of the simulation domain or in two levels with a single VIdx variable providing the unit box index the particle belongs to and a single VReal variable holding the offset from the center of the unit box This functions changes the format from the two level format to the single level format the updated position is stored in pos 3 3 3 changePosFormatlLvTo2Lv void changePosFormatlLvTo2Lv const VReal amp pos Vidx amp vidx VReal amp vOffset Particle position can be represented in a single level with a single VReal varaible storing the offset from the low end corner of the simulation domain or in two levels with a single VIdx variable providing the unit box index the particle belongs to and a single VReal variable holding the offset from the center of the unit box This functions changes the format from the single level format to the two level format the updated position is stored in v dx and vOffset 3 4 Macros 3 4 1 CHECK CHECK expression Abort if expression is false if this macro is enabled in compile time Section 6 2 1 3 4 2 OUTPUT OUTPUT m
84. pecies elemidx should be equal to or larger than 0 and smaller than the number of molecular species in the simulation instance e std string name Set the name of the molecular species e bce_tyep_e aa_bcType 3 2 Boundary condition types for the x aa_bcType 0 y aa_ bcType 1 and z aa_bcType 2 directions and the low aa_bcType 0 and high aa_ bcType 1 sides Irrelevant if periodic boundary condition is applied e REAL aa_bcVal 3 2 Set the constant boundary value for the low and high sides in the x y and z directions Relevant only when the PDE boundary type is set to BC TYPE DIRCHLET_CONST or BC_TYPE_NEUMANN_CONST e REAL errorThresholdVal If any of the PDE variables becomes smaller than errorThresh oldVal Biocellion prints an error message and aborts the program execution errorThresh oldVal should be equal to or smaller than 0 0 Set errorThresholdVal to REAL_MAX de fined in SBIOCELLION_ROOT include base_type h x 1 0 to turn this check off e REAL warningThresholdVal If any of the PDE variables becomes smaller than warningTh resholdVal but still no variable is smaller than errorThresholdVal Biocellion prints a warn ing message but does not abort the program execution warningThresholdVal should be equal to or smaller than 0 0 and should be equal to or larger than errorThresholdVal Set warning ThresholdVal to REAL MAX defined in SBIOCELLION ROOT include bas e_type h
85. plusplus com doc tutorial classes explains C class class object or instance and access specifiers private and public e http www cplusplus com reference vector vector explains C stan dard template library STL vector e http en wikipedia org wiki Class_variable explains static member vari ables and functions e http en wikipedia org wiki C 2B 2B template explains C template e http en wikipedia org wiki Reference_ C 2B 2B explains the concept of reference in C 67 9 0 Bibliography Colella P DT Graves JN Johnson HS Johansen ND Keen TJ Ligocki DF Martin PW Mc Corquodale D Modiano PO Schwartz TD Sternberg and BV Straalen 2012 Chombo Software Package for AMR Applications Design Document Lawrence Berkeley National Laboratory Daly MJ 2000 Engineering radiation resistant bacteria for environmental biotechnology Cur rent Opinion in Biotechnology 11 3 280 285 Ha SJ JM Galazka SR Kim JH Choi X Yang JH Seo NL Glass JHD Cate and YS Jin 2011 Engineered Saccharomyces cerevisiae capable of simultaneous cellobiose and xylose fermenta tion Proceedings of the National Academy of Sciences 77 16 5822 5825 Hall Stoodley L JW Costerton and P Stoodley 2004 Bacterial biofilms from the Natural environment to infectious diseases Nature Reviews Microbiology 2 95 108 Intel Corporation 2008 Intel Ordinary Differential Equation Solver Lib
86. ption is that cells move divide and die at moderate rates Cell state changes affect the secretion and uptake rates of extra cellular molecules Molecular concentration changes in the extra cellular space drive cell state changes as well Biocellion al lows users to couple these two processes more tightly Users can split a baseline time step into one or more state and grid time steps The cell state update module and the extra cellular space state update module communicate once every state and grid time step A single state and grid time step can be further partitioned to solve partial differential equations PDEs The next section Section 1 4 discusses time steps used to solve PDEs in more detail 1 4 Partial Differential Equation PDE Solvers The current version of Biocellion provides solvers for three different types of PDEs steady state linear reaction diffusion PDEs time dependent linear reaction diffusion PDEs and time dependent reaction diffusion PDEs adopting the splitting scheme Strang 1968 to update the re action part Future releases of Biocellion will support advection reaction diffusion PDEs Biocellion uses CHOMBO Colella et al 2012 as an internal solver for PDEs and uses the Intel ordinary differential equation ODE solver Intel Corporation 2008 to solve the reaction part of the splitting PDE We briefly describe PDE solver features relevant to Biocellion users in this section and refer readers to Colella et al
87. r the elemIdx molecular species This model routine is invoked in ev ery PDE time step and is relevant only for time dependent linear PDEs PDE TYPE REACTION DIFFUSION_TIME_DEPENDENT_LINEAR This model routine is called only when users con figure to call this function Section 2 4 12 vier indexes this unit box ubAgentData holds the entire set of discrete agents in this unit box gridPhi is the concentration of the elemIdx molec ular species at the beginning of a PDE time step v_gridModelReal and v_gridModelInt are the model specific REAL type variables and model specific S32 type variables for this unit box re spectively gridRHS is initially set to the gridRHS value set inside the updatelfGridRHSLinear function for this molecular species Users may change the value based on model specifics 50 4 5 10 updatelfGridRHSTimeDependentSplitting void updatelfGridRHSTimeDependentSplitting const S32 pdeldx const VIidx amp vidx const UBAg entData amp ubAgentData const Vector lt double gt amp v_gridPhi const Vector lt REAL gt amp v_gridMod elReal const Vector lt S32 gt amp v_gridModelInt Vector lt double gt v_gridRHS Set the derivative values for the reaction term of the pdeIdx PDE which is solved using the split ting scheme v dx indexes this unit box ubAgentData holds the entire set of discrete agents in this unit box v_gridPhi holds molecular concentrations for the molecular spec
88. r the face between two adjacent unit boxes model routines set the parameters based on the states of the two adjacent unit boxes and the states of the discrete agents residing in the two boxes Figure 1 5 A model routine updating the state of a discrete agent can access the state of the unit box the agent belongs to red and its 26 neighboring unit boxes white Biocellion allows users to set sub regions of the simulation domain as a partial differential equation PDE buffer the PDE buffer is relevant only when solving PDEs No discrete agent is allowed to reside in the PDE buffer Say we are simulating a yeast colony growing on top of a thick agarose cylinder The agarose cylinder holds nutrients and other molecules secreted by yeast cells The agarose cylinder region needs to be included in simulation to track molecular concentration changes in the extra cellular space with high accuracy However if no cells reside in the agarose cylinder region maintaining data structures for all three computational modules at the resolution of the interface grid spacing can waste a significant amount of computing and memory Users can set sub regions of the simulation domain as a PDE buffer and adopt a coarser grid spacing for the buffer Biocellion supports adaptive mesh refinement AMR Section 1 4 and the PDE buffer grid spacing will be matched to the coarsest AMR level grid spacing Figure 1 6 provides an example
89. rary Reference Manual Jiao Y and S Torquato 2011 Emergent Behaviors from a Cellular Automaton Model for Invasive Tumor Growth in Heterogeneous Microenvironments PLoS Computational Biology 7 12 e1002314 Kang S S Kahan and B Momeni TBP Methods in Molecular Biology Ch Simulating Microbial Community Patterning using Biocellion Springer Kim Y MA Stolarska and HG Othmer 2007 A hybrid model for tumor spheroid growth in vitro I theoretical development and early results Mathematical Models and Methods in Applied Sciences 17 1 1773 1798 Majors PD JS McLean and JC Scholten 2008 NMR bioreactor development for live in situ microbial functional analysis Journal of Magnetic Resonance 192 1 159 166 Melnicki MR GE Pinchuk EA Hill LA Kucek SM Stolyar JK Fredrickson AE Konopka and AS Beliaev 2013 Feedback controlled LED photobioreactor for photophysiological studies of cyanobacteria Bioresource Technology 134 127 133 Momeni B KA Brileya MW Fields and W Shou 2013 Strong inter population cooperation leads to partner intermixing in microbial communities eLife p 00230 Newman TJ 2005 Modeling Multicellular Systems Using Subcellular Elements Mathemati cal Biosciences and Engineering 2 3 61 1 622 Resat H V Bailey LA McCue and A Konopka 2011 Modeling Microbial Dynamics in Het erogeneous Environments Growth on Soil Carbon Sources Microbial Ecology Rubinstein MR
90. re used to provide model specifics see Section 4 for the user interface to provide model specifics e domain bdry type e Used to set the simulation domain boundary types in the x y and z directions DOMAIN BDRY TYPE NONPERIODIC HARD WALL Non periodic boundary condi tion is assumed In solving PDEs users can apply either Dirichlet or Neumann bound ary conditions for the low and high sides Discrete agents cannot pass the domain boundary DOMAIN BDRY TYPE PERIODIC Periodic boundary condition is assumed for PDEs Discrete agents passing the domain boundary will appear at the other end of the simu lation domain e pde buffer bdry type e Used to set the boundary type between an interface grid region and a PDE buffer region This boundary type is only related to discrete agent movement and irrelevant to PDE solves PDE BUFFER BDRY TYPE HARD WALL This is the only valid option at this point Discrete agents cannot penetrate into a PDE buffer region PDE buffer regions are automatically marked as uninhabitable e pde type e Used to set PDE types PDE TYPE REACTION DIFFUSION STEADY STATE LINEAR Used to find a quasi steady state solution of a reaction diffusion linear PDE This option assumes that the PDE reaches steady state within a single state and grid time step PDE TYPE REACTION D
91. rovides interface functions to access the variables The entire set of discrete agent state variables are stored in an internal linear array Biocellion does not allow direct access to the array which can be error prone Public non static member functions e agentType_t getType void const Return the type of the discrete agent void setType const agentType_t type Set the type of the discrete agent to type This function also resizes and resets the internal linear array Users need to properly initialize discrete agent state variables after this function is invoked e REAL getRadius void const Return the value of the radius variable a default discrete agent state variable e void setRadius REAL radius Set the value of the radius variable a default discrete agent state variable to radius e void getODEVal const S32 odeNetldx const S32 varldx const Return the value of the varldx variable in the odeNetIdx ordinary differential equation ODE system A single discrete agent can have more than one ODE systems 23 e Vector lt REAL gt getODEValArray const S32 odeNetldx const Return a Vector array hold ing the entire set of variable values in the odeNetIdx ODE system e void setODEVal const S32 odeNetldx const S32 varldx REAL val Set the value of the varldx variable in the odeNetIdx ordinary differential equation ODE system to val e void setODEValArray const S32 odeNetldx const Vector lt REAL gt
92. rtVIdx const VIdx amp pdeB ufferBoxSize const REAL gridPhi REAL amp gridRHS Update the source term for the elemidx molecular species This model routine is relevant only for linear PDEs PDE TYPE REACTION DIFFUSION STEADY STATE LINEAR and PDE TYPE REACTION DIFFUSION TIME DEPENDENT LINEAR startVIdx points the low end cor ner of the PDE buffer unit box assuming the interface grid spacing pdeBufferBoxSize is the size of a PDE buffer unit box assuming the interface grid spacing gridPhi holds the molecular concen tration of the elemIdx molecular species for this PDE buffer grid unit box This model routine updates gridRHS 4 5 20 adjustPDEBufferRHSTimeDependentLinear void adjustPDEBufferRHSTime DependentLinear const S32 elemldx const VIdx amp startVIdx con st Vidx amp pdeBufferBoxSize const REAL gridPhi REAL amp gridRHS 53 Adjust the source term for the elemIdx molecular species This model routine is invoked in ev ery PDE time step and is relevant only for time dependent linear PDEs PDE TYPE REACTION DIFFUSION_TIME_DEPENDENT_LINEAR This model routine is called only when users con figure to call this function Section 2 4 12 startVIdx points the low end corner of the PDE buffer unit box assuming the interface grid spacing pdeBufferBoxSize is the size of a PDE buffer unit box assuming the interface grid spacing gridPhi is the concentration of elemI
93. s the entire set of discrete agents in this unit box v_gridPhi v_gridModelReal and v_gridModelInt are the molecular concentrations model specific REAL type variables and model specific S32 type variables for the unit box in the interface grid partition respectively This model routine updates gridBeta 4 5 8 updatelfGridRHSLinear void updatelfGridRHSLinear const S32 elemldx const Vidx amp vid const UBAgentData amp ubAge ntData const Vector lt REAL gt v_gridPhi const Vector lt REAL gt v_gridModelReal const Vect or lt S32 gt 4 v_gridModellnt REAL amp gridRHS Update the source term for the elemIdx molecular species This model routine is relevant only for linear PDEs PDE TYPE REACTION DIFFUSION STEADY STATE LINEAR and PDE_ TYPE REACTION DIFFUSION TIME DEPENDENT LINEAR vidx indexes this unit box ubAgent Data holds the entire set of discrete agents in this unit box v gridPhi v gridModelReal and v_gridModelInt are molecular concentrations model specific REAL type variables and model specific S32 type variables for this unit box respectively This model routine updates gridRHS 4 5 9 adjustIfGridRHSTimeDependentLinear void adjustIfGridRHSTimeDependentLinear const S32 elemldx const VIdx amp vide const REAL e ridPhi const Vector lt REAL gt v gridModelReal const Vector lt S32 gt amp v_gridModelInt REAL amp gridRHS Adjust the source term fo
94. single variable can be updated by multiple model rou tine calls This requires a mechanism to coordinate the multiple updates often referred as data synchronization sync_method_e is used to set the coordination method when a single variable is updated in multiple function calls SYNC_METHOD_OVERWRITE If this option is selected only one of the multiple up 13 dates is reflected e g if three different model routines set the value of variable a to 3 5 and 9 respectively a is set to either 3 5 or 9 with equal probability SYNC METHOD DELTA If this option is selected the differences from the initial value are summed e g if three different model routines set the value of variable a to 3 5 and 9 respectively then a is set to 3 5 9 assuming that the initial value is 0 SYNC_METHOD_TRANSACTIONAL Reserved for future use e if grid var type e Used to set the type of a variable associated with a unit box in the interface grid IF_GRID_VAR_TYPE_PHI Set the type to a REAL variable storing molecular con centration which is mainly updated by solving PDEs Biocellion updates variables of this type both when solving PDEs and when inside the model routine updating grid state variables based on model specific rules Model specific REAL and S32 variables are updated only inside the model routine updating grid state variables IF_GRID_VAR_TYPE_MODEL_REAL Set the type to
95. st image left is reproduced from http en wikipedia org wiki Yeast Updating individual cell states is largely a users task at this point Biocellion provides a solver for deterministic ordinary differential equation ODE systems but there are many other approaches to modeling cell regulatory networks We plan to add solvers for widely used regulatory network modeling approaches e g boolean network stochastic ODE systems and hybrids For a cell mapped to multiple discrete agents we can consider each discrete agent as a cell compartment or a subcellular organelle Another path to extend Biocellion is to provide solvers to update the state of subcellular discrete agents or to support modeling of molecular transport across subcellular compartments Computing the flow rate of liquid in biological systems with time varying structures e g water flow in soil aggregate or blood flow in human tissue is another important computational task in building high fidelity models In computing blood flow rate our current interest is mainly in tracking long term average flow rate changes in capillaries which affect the delivery rate of nutrients and vessel stabilization and regression to model tissue or vascular tumor growth rather than modeling short term dynamics e g flow changes within a single heart beat cycle of blood flow in arteries which is more relevant in studying cardiovascular disease 2 0 Biocellion Data Types 2 1 Basic Data Types
96. state variables associated with an interface grid unit box v dx indexes the unit box ubAgentData holds the discrete agents located in this unit box This model routine updates the molecular concentrations v_gridPhi the model specific REAL type attributes v_gridModelReal 47 and the model specific S32 type attributes v_gridModelInt for this unit box 4 5 2 updatelfGrid Var void updatelfGridVar const BOOL pre const S32 round const VIdx amp vIdx const NbrBox lt const UBA gentData gt 4 ubAgentDataPtrNbrBox Vector lt NbrBox lt REAL gt gt amp v_gridPhiNbrBox V ector lt NbrBox lt REAL gt gt amp v_gridModelRealNbrBox Vector lt NbrBox lt 832 gt gt amp v_gridModel IntNbrBox Update the state variables associated with an interface grid unit box and its 26 neighboring boxes It is an error to access or update values associated with a unit box location outside the simulation domain or a unit box location in a PDE buffer region pre indicates whether this routine is called at the beginning of a state and grid time step if pre is true or at the end of a state and grid time step if pre is false round indicates the round in which this model routine is called users are al lowed update extra cellular space state variables in multiple rounds Section 2 4 6 vidx indexes this unit box ubAgentDataPtrNbrBox holds pointers for UBAgentData class instances for this unit box and its 26 neighboring unit boxes This model routine upda
97. static member variables 17 REAL durationBaselineTimeStep Set the size of baseline time step e S32 numStateAndGridTimeStepsPerBaseline The state and grid time step size is set to baseline time step size divided by numStateAndGridTimeStepsPerBaseline 2 4 6 OptModelRoutineCallInfo Biocellion users can update the interface grid state variables based on model specific rules Bio cellion optionally invokes the model routine updating interface gird state variables at the beginning and at the end of every state and grid time step This class instance controls the optional model routine invocation Public non static member variables e S32 numUpdatelfGridVarPreStateAndGridStepRounds Call the model routine once for every unit box in the interface grid for numUpdatelfGridVarPreStateAndGridStepRounds rounds at the beginning of every state and grid time step This variable should be zero or a positive integer e S32 numUpdatelfGridVarPostStateAndGridStepRounds Call the model routine once for every unit box in the interface grid for numUpdatelfGridVarPostStateAndGridStepRounds rounds at the end of every state and grid time step This variable should be zero or a posi tive integer 2 4 7 SpAgentInfo This class instance holds the model specifics about a discrete agent type Public non static member variables e REAL dMax Set the maximum direct physico mechanical interaction distance The value of this variable should be equal to or smal
98. stems automating for example the parallel solution of partial differential equations and propagation of physical contact forces while also providing developers easy interfaces to express model specifics such as growth rates and ad hesion properties Biocellion provides model developers with flexibility in expressing model specifics while relieving them of the burdens of parallel computation and optimization Today modelers express these model specifics as program functions that Biocellion invokes when it runs This provides modelers with the flexibility of a general purpose programming language when expressing model specifics without significant loss in computational efficiency That is model developers express model specifics as sequential code C code using only basic C features Section 8 explains C concepts relevant to Biocellion The Biocellion core frame work links to this C model library at runtime invoking functions as needed to integrate the model specifics Realizing that most computational biologists are not C experts we en vision adding a productivity layer on top of the current version of Biocellion once mathematical modeling of biological systems matures sufficiently Other possibilities are enabling modelers to use other languages eg Python with which they may be more comfortable at some cost to performance but no loss in scalability 1 1 Simulation Components Biocellion simulates biological systems with a large n
99. tep load balance once per every 100 baseline time steps regrid AMR hierarchy once per every 100 baseline time steps and generate checkpoint data every 600 baseline time steps The load balancing interval should be a positive integer multiple of the regirdding interval Setting the interval to 0 turns off summary load balancing regridding or checkpoint data generation Default values for summary output and checkpoint data generation are 0 and load balancing and regirdding are 100 e amr Set the refinement ratio between two consecutive AMR levels and the desired fill ratio for AMR boxes The fill ratio is the number of user tags in an AMR box divided by the number of unit boxes in the AMR box If fill ratio is high Biocellion tends to create a large number of tiny AMR boxes which lowers the efficiency of computation If fill ratio is set to a low value Biocellion tends to create fewer boxes but a larger percentage of the simulation domain is covered by fine grids The default refinement ratio is 2 and the refinement ratio should be either 2 or 4 The default fill ratio is 0 5 and 0 0 lt fill ratio lt 1 0 6 4 Execution on Desktop PCs and Workstations Biocellion usees Intel Thread Building Block version 4 2 or later for multi threading and users need to set the LD LIBRARY PATH Linux environment variable to point the TBB library direc tory in TBB 4 2 this is STBB_ROOT 1lib intel64 gcc4 1 Biocellion executables are located under SBIO
100. tes the molecular concen trations v_gridPhiNbrBox the model specific REAL type attributes v_gridModelRealNbrBox and the model specific S32 type attributes v gridModelIntNbrBox for this unit box and its 26 neighboring unit boxes v_gridPhiNbrBox v_gridModelRealNbrBox and v_gridModelIntNbrBox initially holds the current values associated with this box and its 26 neighboring boxes and users may update a subset of the values based on their model specifics 4 5 3 updatelfGridKappa void updatelfGridKappa const Vidx amp vIdx const UBAgentData amp ubAgentData const Vector lt REAL gt amp v_gridPhi const Vector lt REAL gt amp v_gridModelReal const Vector lt S32 gt amp v_gridMo delInt REAL amp gridKappa Set the ratio of the available extra cellular volume in computing diffusion flux molecular concen tration in Biocellion is the total amount of a molecular species in a unit box divided by the volume of the unit box However in computing diffusion the direction and rate of diffusion is determined by the molecular concentration in the free space If one unit box is heavily occupied by cells and another unit box has no cell then the net amount of diffusion to the empty unit box can be positive even when the molecular concentration ignoring cell volume exclusion in the empty unit box is higher This model routine sets gridKappa to consider cell volume exclusion in computing diffu sion flux Say there are two unit boxes box0 and bo
101. the absolute value of the solution is larger than threshold relative error is controlled Otherwise absolute error threshold x epsilon is controlled 2 4 12 PDEInfo This class instance holds the information required to solve a PDE Public non static member variables e pde_type_e pdeType Set the type of the PDE e S32 numLevels Set the number of adaptive mesh refinement AMR levels e S32 numTimeSteps Set the number of PDE time steps per state and grid time step e BOOL callAdjustRHSTimeDependentLinear If set to true Biocellion invokes the model routines adjustIfGridRHS Time DependentLinear and adjustPDEBufferRHSTimeDependent Linear at the beginning of every PDE time step to fine tune the PDE reaction term If set to false Biocellion updates the PDE reaction part only once per state and grid time step Relevant only when pdeType is set to PDE TYPE REACTION DIFFUSION TIME DEPENDENT LINEAR e MGSolvelnfo mgSolvelnfo Set the PDE multigrid solver options See Section 2 4 10 for additional details e AdvectionInfo advectionInfo Reserved for future use Irrelevant at this point e SplittingInfo splittingInfo Set the information required to solve the PDE using the splitting scheme Relevant only when pdeType is set to PDE TYPE REACTION DIFFUSION TIME DEPENDENT SPLITTING 21 e Vector lt GridPhilnfo gt v_gridPhilnfo Provide the inf
102. the discrete agent is a sphere vOffset is the offset from the center of the unit box containing the discrete agent and points the center of the sphere and radius is the radius of the sphere The unit box containing the discrete agent and its 26 neighboring boxes are considered the radius of a discrete agent cannot exceed one half of the interface grid spacing and the sphere can overlap with only these 3 x 3 x 3 unit boxes The first function finds the ratios for the entire set of 3 x 3 x 3 unit boxes aaa ratio contains the solution and the second function returns the ratio for only one of the 3 x 3 x 3 unit boxes indexed by ubVIdxOffset This function partitions the minimum bounding box of the sphere to 2 Level x qmaxLevel x gmaxLevel sub boxes We consider that a sub box belongs to a unit box if the sub box center is located inside the sphere and the unit box If the center is located exactly at the surface of the sphere we consider the center is located inside the sphere If the center is located exactly at the boundary of two unit boxes we consider that the center is located at the unit box with a larger index These functions compute the ratio by dividing the 33 number of sub boxes belong to a unit box by the number of sub boxes belong to the 3 x 3 x 3 boxes maxLevel should be a non negative integer equal to or smaller than 7 A larger value increases both accuracy and execution time vOffset cannot point the location outside the unit b
103. to respectively dir is a unit vector from the position of spAgent to the position of spAgent0 dist is the distance between the positions of the two discrete agents This model routine updates the temporary extra mechanical interaction states for spAgent0 and spAgent which are extraMechIn trctData0 and extraMechIntretDatal respectively Temporary mechanical interaction state vari ables are reset to zero at the beginning of every baseline time step Users also form and break a junction between a discrete agent pair in this model routine If the two discrete agents are linked and link is set to true the junction ends in the spA gen and spAgentl sides are replaced with end and endl respectively isLinked Section 2 4 19 returns whether two discrete agents are linked or not If spAgent0 and spAgentl are linked and unlink is set to true the Biocellion framework breaks the junction If there is no junction between spAgent0 and spAgentl and link is set to true a new junction is formed with junction ends end and end If there is no junction and unlink is set to true or both link and unlink is set to false no change occurs It is an error to set both link and unlink to true 4 5 State Changes in the Extra cellular Space 4 5 1 initIfGridVar void initlfGridVar const VIdx amp vIdx const UBAgentData amp ubAgentData Vector lt REAL gt v_g ridPhi Vector lt REAL gt amp v_gridModelReal Vector lt 832 gt amp v_gridModellnt Initialize the
104. to true if this discrete agent is subject to disappear If neither divide nor disappear is set to true adjustSpA gent Section 4 3 6 is called It is an error to set both divide and disappear to true updateSpA gent BirthDeath is called after spA gentSecretionBySpAgent within a single baseline time step 4 3 6 adjustSpAgent void adjustSpA gent const VIdx amp vidx const AgentJunctionInfo amp junctionInfo const VReal amp vO Deet const AgentMechIntrctData amp mechIntrctData const Vector lt NbrBox lt REAL gt gt amp v_gridP hiNbrBox const Vector lt NbrBox lt REAL gt gt amp v_gridModelRealNbrBox const Vector lt NbrBox lt 32 gt gt amp v gridModelIntNbrBox SpAgentState amp state VReal amp disp Adjust the state and position of this discrete agent v dx indexes the unit box this discrete agent is located in junctionInfo stores the list of junctions between this discrete agent and its neigh boring discrete agents vOffset is the offset from the center of the unit box this discrete agent be longs to mechIntrctData stores the temporary mechanical interaction state of this discrete agent v_gridPhiNbrBox v_gridModelRealNbrBox and v_gridModelIntNbrBox provide the information molecular concentrations model specific REAL type attributes and model specific S32 type at tributes respectively about the extra cellular space for this box and the 26 neighboring boxes This model routine updates the state of this discrete agent
105. tra output variables in addition to the default radius and color variables for each particle representing a discrete agent This value is set inside the simulation configuration file see Section 6 3 3 1 11 getModelParam std string getModelParam void 32 Return the model parameter string set inside the simulation configuration file see Section 6 3 3 1 12 getAMRRatio S32 getAMRRatio void Return the refinement ratio between two consecutive adaptive mesh refinement AMR levels The refinement ratio is set inside the simulation configuration file see Section 6 3 3 2 Random Number Generator 3 2 1 getModelRand REAL Util getModelRand const S32 idx Return a random number generated using the idx random number generator Biocellion allows users to assign a set of random number generators with different probability density functions 3 3 Extra Support Functions 3 3 1 computeSphereUBVolOvlpRatio void computeSphereUBVolOvlpRatio const S32 maxLevel const VReal amp vOffset const REAL ra dius REAL aaa ratio 3 3 3 REAL computesSphereUBVolOvlpRatio const S32 maxLevel const VReal amp vOffset const REAL radius const Vidx amp ubVIdxOffset A discrete agent with a non zero volume may span multiple unit boxes Modelers may want to find the ratios of the overlapping volumes between the agent and different unit boxes over the entire volume of the agent These functions find an approximate solution assuming that
106. umber of cells capturing individual cell behaviors interactions between cells and interactions of cells with their biotic and abiotic environment The current version of Biocellion maps each cell to a discrete agent In future versions Biocellion will allow users to map a cell to multiple discrete agents as illustrated in Newman 2005 Sander sius et al 2011 Discrete agents can also represent other non cellular biological elements e g Xavier et al represent extra cellular matrix using discrete agents in Xavier et al 2005 Each discrete agent has associated variables representing its position and state Biocellion users set the number of discrete agent types and the number of discrete agent state variables or attributes for each discrete agent type Each agent can be placed anywhere in the simulation domain unless restricted by the mathematical model to within the limits of floating point precision Individ ual agents move divide disappear and secrete new discrete agents Biocellion considers each discrete agent as a sphere by default and each discrete agent has radius as a state variable Op tionally attributes may be added to a discrete agent to define a non spherical shape Figure 1 1 provides an example Four attributes three for direction and one for height are added to treat a discrete agent as a cylinder The Biocellion core framework does not use radius or other discrete agent state variables internally to update the state
107. updated data 4 2 15 init void init void Users may perform model specific initialization inside this function 4 2 16 term void term void 42 Users may perform model specific clean up inside this function 4 2 17 setPDEBuffer void setPDEBuffer const VIdx amp startVldx const VIdx amp regionSize BOOL amp isPDEBuffer Mark whether this partition is an interface grid partition set isPDEBuffer to false or a PDE buffer partition set isPDEBuffer to true startVldx points the low end corner of this partition and regionSize 1s the size of this partition assuming the interface grid spacing 4 2 18 setHabitable void setHabitable const VIdx amp vidx BOOL amp isHabitable Mark whether this unit box indexed by vidx is habitable set isHabitable to true or not set isHabitable to false If this unit box is marked as uninhabitable discrete agents cannot penetrate into this unit box 4 3 Individual Agent Behavior 4 3 1 addSpAgents void addSpA gents const BOOL init const Vidx amp startVIdx const VIdx amp regionSize const IfGrid BoxData lt BOOL gt amp ifGridHabitableBoxData Vector lt VIdx gt 4 v_spAgentVIdx Vector lt SpAge ntState gt amp v_spAgentState Vector lt VReal gt amp v_spAgentOffset Add discrete agents to this partition at the beginning of a simulation if init is true or at the beginning of a baseline time step if init is false startVldx points the unit box at the low end corner of this p
108. variables for the unit box spAgent belongs to respectively dir isa unit vector from the position of spAgent to the position of spAgent0 dist is the distance between the positions of the two discrete agents This model routine updates force 4 4 3 computeExtraMechIntrctSpAgent void computeExtraMechIntrctSpAgent const Vidx amp vIdx0 const SpAgent amp spAgent0 const Vect or lt REAL gt amp v_gridPhi0 const Vector lt REAL gt v_gridModelReal0 const Vector lt S32 gt amp v_g ridModellnt0 const VIdx amp vIdx1 const SpAgent amp spAgentl const Vector lt REAL gt v_gridPhi I const Vector lt REAL gt amp v_gridModelReall const Vector lt S32 gt 4 v_gridModelInt1 const VRe al amp dir const REAL amp dist ExtraMechIntrctData amp extraMechIntrctData0 ExtraMechIntrctData amp extraMechIntrctDatal BOOL amp link JunctionEnd amp end0 JunctionEnd amp endl BOOL amp unlin k Evaluate extra mechanical interactions between a discrete agent pair v dx0 indexes the unit box spAgent0 is located in v_gridPhi0 v_gridModelReal0 and v_gridModelIntO are molecular con centrations model specific REAL type variables and model specific S32 type variables for the unit box spAgent0 belongs to respectively v dx indexes the unit box spAgentl is located in v_gridPhil v_gridModelReall and v_gridModelInt1 are molecular concentrations model specific REAL type variables and model specific S32 type variables for the unit box spA gent belongs
109. x1 sharing a face and dy are molecular concentrations for box0 and box respectively ko and x are gridKappa values set for box0 and box1 respectively is the diffusion coefficient set for the common face Then the diffusion flux from box0 to box becomes B x do Y x i vidx indexes this unit box ubAgentData holds the entire set of discrete agents in this unit box v_gridPhi v_gridModelReal and v_gridModelInt are molecular concentrations model specific REAL type variables and model specific 32 type variables for this unit box respectively 0 0 lt gridKappa lt 1 0 Set gridKappa to 1 0 to ignore cell volume exclusion gridKappa is assumed to be 1 0 outside the simulation domain 48 4 54 updatelfGridAlpha void updatelfGridAlpha const S32 elemldx const VIdx amp vIdx const UBAgentData amp ubAgentDa ta const Vector lt REAL gt amp v_gridPhi const Vector lt REAL gt amp v_gridModelReal const Vector lt S32 gt amp v_gridModellnt REAL amp gridAlpha Set the decay rate for the elemIdx molecular species vidx indexes this unit box ubAgent Data holds the entire set of discrete agents in this unit box v_gridPhi v_gridModelReal and v_gridModellnt are molecular concentrations model specific REAL type variables and model specific S32 type variables for this unit box respectively This model routine updates gridAlpha 1 0 lt gridAlpha lt 0 0 Set gridAlpha to 0 0 to ignore decay 4 5 5 updatelfGridBeta

Biocellion 1.1 User Manual

Contents

Download Pdf Manuals

Related Search

Related Contents