Digitize your Biology! Modeling multicellular systems through interpretable cell behavior

Cells are fundamental units of life, constantly interacting and evolving as dynamical systems. While recent spatial multi-omics can quantitate individual cells’ characteristics and regulatory programs, forecasting their evolution ultimately requires mathematical modeling. We develop a conceptual framework—a cell behavior hypothesis grammar—that uses natural language statements (cell rules) to create mathematical models. This allows us to systematically integrate biological knowledge and multi-omics data to make them computable. We can then perform virtual “thought experiments” that challenge and extend our understanding of multicellular systems, and ultimately generate new testable hypotheses. In this paper, we motivate and describe the grammar, provide a reference implementation, and demonstrate its potential through a series of examples in tumor biology and immunotherapy. Altogether, this approach provides a bridge between biological, clinical, and systems biology researchers for mathematical modeling of biological systems at scale, allowing the community to extrapolate from single-cell characterization to emergent multicellular behavior.


Cycling
As introduced in prior work 2,13,14 , a cell cycle model is a sequence of phases (with indices  = 0,1, … , ) and transition rates (or exit rates) { !} connecting phase  to the next phase  + 1.By convention, exit from the first phase (with index 0) is called cycle entry, and cells divide upon exiting last phase (with index ) to return to the 0 th phase.In the reference implementation, cell transitions between phases are probabilistic, based upon the transition rate; if a cell is in phase  at time , then the probability of advancing to phase  + 1 between the current time  and a future time  +  is given by (exit phase ) =  !Δ. (1) Note that the mean time spent in the  th phase is For example, a five-phase model ("flow cytometry separated") consists of phases G0, G1, S, G2, and M, cycle entry is a transition from G0 (Phase 0) to G1 (Phase 1), and division occurs at the end of M (Phase 4).In the reference implementation 13 , at division daughter cells are assigned half the volume of the parent cell and positioned randomly such that (1) they both fit within the parent cell's volume, and (2) the daughter cells have the same center of volume as their parent cell.In the current grammar, the dictionary of cell cycling behaviors and their corresponding parameters are in Table 1.Future versions of the grammar may include finer grained control on placement of daughter cells after division.

Death
In the initial grammar, we support apoptosis (as a primary form of non-immunogenic cell death) and necrosis (as a key form of immunogenic cell death).We note that at a cellular and multicellular perspective, death is not merely a discrete event, but rather entry into a cascade of processes [15][16][17][18][19][20] ; as noted in prior modeling and analysis 13,15,21 , every model of death has a nonzero, finite duration on the order of hours (for apoptosis) to days or weeks (for necrosis).In the reference implementation 13 , live cells have an apoptotic death rate  $ , such that the probability that a cell becomes apoptotic between the current time  and a future time  + Δ is given by (cell becomes apoptotic) =  $ Δ. (2) Once the cell becomes apoptotic, it shrinks until being removed from the simulation.
Similarly, live cells have a necrotic death rate  % , such that the probability that a cell becomes necrotic between the current time  and a future time  + Δ is given by (cell becomes necrotic) =  % Δ (3)   dictionary symbol synonyms controllable phenotype parameter exit from cycle phase 0 cycle entry  " exit from cycle phase 1  # exit from cycle phase 2  $ exit from cycle phase 3  % exit from cycle phase 4   & exit from cycle phase 5  ' TABLE 1. Cell cycle terms in the behavior dictionary, along with corresponding controllable parameters in the PhysiCell reference implementation.
In the reference necrosis model, cells initially swell until reaching a bursting volume (by default, 200% of the cell's volume at the start of necrosis), and then gradually shrink over the course of days.
Full details on the reference parameter values for apoptosis and necrosis (including rates of cell volume change and durations) are found in prior work 13,21 .In the current grammar, the dictionary of behaviors and their corresponding parameters are in Table 2.

Secretion, uptake, and generalized chemical export
Cells secrete and uptake (consume) chemical factors that diffuse in the extracellular microenvironment, with a fixed mathematical form first introduced in BioFVM 22 and later incorporated in PhysiCell 13 .For any diffusible substrate , BioFVM solves the reaction-diffusion equation in Equation ( 4): Here,  ! is the  th cell's position, () is the Dirac delta function that mathematically focuses the cellbased source/sink at its center,  ! is its secretion rate of the substrate (with dimensions 1/time),  !* is the "target value" of its secretion (i.e., secretion slows as  approaches this value), and  ! is its uptake or consumption rate (also with dimensions 1/time).When these fixed secretion and uptake forms are insufficient to match experiments or advanced modeling forms, we provide a generic net export term  !(with dimensions substrate/cell/time).Because the microenvironment can have multiple diffusing substrates, the grammar can automatically expand to reference secretion, uptake, and export parameters for each one, therefore defining a family of symbols as summarized in Table 3.

(Biased) migration and chemotaxis
As previously introduced in PhysiCell 13 , cell migration (or motility) is represented as a biased random walk: the cell's migration direction  7., is a combination of a random (unit vector) component  and a non-random directed component (a migration bias direction)  8-2( , which is normalized to be a unit vector as in Equation ( 5):  3. Secretion, uptake, and export terms in the behavior dictionary, along with corresponding controllable parameters in the PhysiCell reference implementation.X denotes the name of any diffusible substrate (e.g., oxygen).
Here, 0 ≤  ≤ 1 is the cell's migration bias: if  = 1, then migration occurs completely along the bias direction  8-2( , while  = 0 corresponds to purely Brownian motion.We then obtain the migration velocity  7., by multiplying the direction by the migration speed (s): In our formulation, the cell's migration speed (and hence migration velocity) can change dynamically based upon any further modeling rules, even if the migration direction does not.Cell migration has a persistence time  1)+(-(, : between the current time  and a future time  + Δ, the probably of choosing a new migration direction (by re-evaluating Equation ( 6)) is given by Note that this gives a mean time of  1)+(-(, between direction changes.
Chemotaxis is modeled by setting the migration bias direction along the gradient of one or more chemical substrates in the microenvironment.If available chemical substrates are  9 ,  " , … ,  : , then we set where  ! is a weighting (chemotactic response) for each chemical gradient:  !< 0 for migration against the gradient (along −∇ ! ),  !> 0 for migration along the gradient, and  != 0 if the gradient makes no contribution.Because the microenvironment can have multiple diffusing substrates, the grammar can automatically expand to a family of chemotactic response terms.See Table 4.

Cell-cell adhesion
Two models of cell-cell adhesion are supported: a looser cell-cell adhesion using potential functions (as in prior mathematical models [23][24][25][26][27][28] ), and more persistent elastic springs 6,29,30 .In the reference PhysiCell implementation for potential-based adhesions, if cells  and  are adhered and within interaction distance, then the contribution to the velocity of cell  is given by where  ! is cell 's adhesive strength,  !; is the (relative) adhesive affinity of cell  to cell  (which for example could be modeled based upon adhesive receptor expressions of cells  and ),  ! is the position of cell , and  $,! is cell i's relative maximum adhesion distance (the largest distance it can extend for adhesion, as a multiple of its radius), and  ! is the cell's radius 13 .Notice that this adhesion only depends dictionary symbol synonyms controllable phenotype parameter migration speed  migration bias  migration persistence time  ,-./0/1 chemotactic response to X chemotactic sensitivity to X  * TABLE 4. Migration terms in the behavior dictionary, along with corresponding controllable parameters in the PhysiCell reference implementation.X denotes the name of any diffusible substrate (e.g., oxygen).
upon the relative distance between cells, and not their prior history.If the cells are identical, the coefficient reduces to  !!  ! .See Table 5.
If more sophisticated cell-cell adhesion is required (e.g., to model that cell adhesions may form more readily than they break), we also support an elastic cell-cell adhesion model.Cell  forms an elastic adhesion to cell  between the current time  and a future time  + Δ at an attachment rate  $,!; (and thus probability  $,!; Δ), provided that: • The distance ! !−  " ! between the cells does not exceed the maximum adhesion interaction distance, and • Neither cell  nor cell  has exceeded their individual maximum number of cell adhesions ( #," and  #,! ).
We calculate this cell-cell attachment rate  $,!; as where  $,! is cell 's rate of forming elastic attachments and  !; is the adhesive affinity of cell  to cell .(Notice that cell  can form attachments to cell , and cell  can form attachments to cell  independently with independent rates.)Analogously, cell  can detach from cell  at rate  ?,! , so that during any time interval  to  + Δ, the cell  can detach from cell  with probability  ?,! Δ.As previously modeled [ref], when the cells are elastically adhered, the adhesion contributes to cell  velocity via: where  ! is cell i's elastic constant.Notice that if cells i and j are identical, then the coefficient simplifies to  !;  ! .The full parameter list for basic and elastic cell-cell adhesion (and cell repulsion) is found in Table 5.

Resistance to deformation and movement
PhysiCell 13 , as many agent-based simulation frameworks [23][24][25][26][27] , uses cell "repulsion" to model resistance to deformation, compression, and cell overlap.In the reference PhysiCell implementation for potential-based repulsion, if cells  and  are in contact, then the contribution to the velocity of cell  is given by  5. Cell-cell adhesion and repulsion terms in the behavior dictionary, along with corresponding controllable parameters in the PhysiCell reference implementation.Here, X denotes any cell type.
where  ! is cell i's radius and  ! is its cell-cell repulsion coefficient as defined in prior work [ref].In cases where a cell should be regarded as a rigid object or obstacle (e.g., if rigidly adhered to an underlying matrix), it can be flagged as unmovable using a "movable" parameter m.If m is false ( = 0), then the cell can impart adhesive and repulsive forces on other cells, but they cannot impose reciprocal forces that contribute to the cell's movement.The full parameter list for cell-cell adhesion and repulsion is found in Table 5.

Transformation (changing type)
Many biological systems require cell type transformations from one type to another, particularly via differentiation, transdifferentiation, and mutation.Mathematically, these all can be represented as a transformation (type change) rate  @,!; from type i to type j, with the probability of a type change between the current time  and a future time  + Δ given by  @,!; Δ.(13)   In the simplest reference implementation, when a cell transforms, it retains its volume and position, and overwrites all other phenotype parameters from the new cell type.These parameters are summarized in Table 6.

Fusion
As a simple reference model of cell fusion, if cell i is in contact with cell j, the cells can fuse between the current time  and a future time  + Δ with probability where  A,!; is the fusion rate for the cell i to the type of cell j.In our reference model, when cells i and j fuse, the newly fused cell: • is placed at the center of volume of the pre-fused cells i and j • combines the volumes of the cells i and j • combines the number of nuclei in cells i and j (if tracked) • combines all internalized substrates in cells i and j (if tracked) At the present time, there is no community consensus on what phenotype parameters should be acquired by the newly fused cell; future versions of the reference model may take a (volume-weighted) average of the pre-fused cells' phenotypes.Moreover, community consensus is required to determine if the fused cell  6.Key cell-cell interaction terms (focused on transformation, fusion, phagocytosis, and effector attack) in the behavior dictionary, along with corresponding controllable parameters in the PhysiCell reference implementation.Here, X denotes any cell type.
should shrink towards a typical single cell size or retain its larger size in the long term.Behavioral parameters that can be modulated by our grammar are summarized in Table 6.

Phagocytosis (or predation or ingestion)
We include a simple reference model of phagocytosis (or predation or ingestion), based on recent modeling 4,5,8 : if cell i is in contact with (live) cell j, then between the current time  and a future time  + Δ, cell i has a probability of phagocytosing cell  given by:  BC,!; Δ, (15)   where  B,!; is cell i's rate of phagocytosing live cells of type j.Similarly, if cell j is dead, then the probability of phagocytosing the cell in that time interval is  B?,! Δ.
In the reference model, when a cell phagocytoses (or ingests) another, it absorbs all its volume, while retaining its original position.(In the event that a more detailed volume is represented as in prior work 13 , cell j's solid volume is added to cell i's cytoplasmic solid volume, and cell j's fluid volume is added to cell i's overall fluid volume.)If the simulation framework actively regulates the volume of cell i (e.g., to maintain a target volume 13 ), then over time cell i will return to its previous volume, as a simplified model of degradation of the phagocytosed cell materials.These parameters are summarized in Table 6.

Effector attack
Based on recent models 4,5,8,31,32 , we define a reference model of cytotoxic effector attack (e.g., CD8 T inflicting fatal damage on a target cell via perforin and granzymes [33][34][35] ).If cell i is in contact with cell j, then its probability of damaging cell j between time  and a future time  +  is given by  $,!; ⋅  ;! Δ, (16)   where  $,!; is cell i's rate of attacking cells of type j, and  ;! is cell j's immunogenicity to cell i.When cell i attacks cell j, the damage  ; of the cell increases by  D272E) Δ ( FG:GHI is the rate at which cell i causes damage to cell j, taken to be 1 in the reference implementation), and the total attack time logged by cell j is increased by Δ.Note that if multiple cells are attacking cell j simultaneously, then both the damage and total attack time increase more rapidly.We note that in this reference model, effector attack increases damage in the target cell j, but it does not cause death without additional hypotheses relating damage to a death rate.Behavioral parameters that can be modulated by our grammar are summarized in Table 6.

Other terms
Because scientists may require custom biology not yet supported in our grammar or reference implementation, we allow a custom symbol to that can access any custom cell parameters.We also reserve symbols for cell-basement mechanical interactions that currently lack reference implementations.See Table 7.
dictionary symbol synonyms controllable phenotype parameter custom:X custom: X, custom X X cell-BM adhesion None (symbol reserved) cell-BM repulsion None (symbol reserved) Table 7.Other cell terms in the behavior dictionary, along with corresponding controllable parameters in the PhysiCell reference implementation.Here, X is a custom cell variable (parameter).
Some signals were included at the request of the mathematical modeling community; for instance, the time signal can trigger global model events (e.g., exposure to a dose of radiation therapy), and apoptotic and necrotic allow more fine-grained control of dead cells.The damage signal is currently a flexible, generic term that could describe membrane, nuclear, mitochondrial, DNA, or other types of damage as needed by the specific context of the modeling application.The total attack time signal similarly allows "area-under-the-curve" (AUC) models that require total exposure time to attacking cells.The custom family of signals enables greater extensibility by allowing direct access to model-specific cell variables, without the need to formally extend the language and its vocabulary itself.
Over time, new signals can readily be added to the standard vocabulary based on community feedback, or when new signals emerge as widespread in modeling.For example, more specific types of damage (e.g., DNA damage or mitochondrial damage) may appear widely in some biological domains, thus justifying the addition of multiple, more specific damage signals in the vocabulary.
We note that while some signals can be found in existing ontologies (e.g., many diffusible chemical substrates in ChEBI 57,58 , mechanical pressure in the ontology for physics in biology (OPB) 59 or non-specific pressure in PATO 60 , or non-specific cell-cell contact in the cell behavior ontology (CBO) 1 , no single ontology was found to incorporate our broader family of signals in a form that can be concretely and algorithmically generated for a specific simulation or mathematical model system.Here, [signal S] and [behavior B] are defined symbols defined above.This construction allows us to efficiently "bundle" multiple behavioral response statements for a cell type.

Optional parameters 1: base and maximum effect
These parameters are used to specify the maximum effect of the rule.They take the form: Here, [base] is the unperturbed value of the behavioral parameter in the absence of signals, and [value] is its saturated response under maximum signal.The values should be stated with units, with a preference for microns for spatial units and minutes for time units.

Optional statements: effect on dead cells
By default, we assume that the hypothesis rules apply to live cells only, to avoid non-physical behaviors such as nonzero motility or proliferation of dead cells.However, rules can be designated to apply to dead cells in addition to live cells with the optional statement: Rule applies to dead cells.

Examples:
Here are examples of behavioral response statements with optional arguments noted in red.
• Pressure decreases cycle entry towards 0 1/min with a linear response, with minimum threshold 0 and maximum threshold 0.5.• doxorubicin increases apoptosis.
• doxorubicin increases apoptosis towards 0.01 1/min with a Hill response, with half-max 0.1 and Hill power 2.

• Virus increases fusion to tumor cells. Rule applies to dead cells.
As an example of bundling multiple statements for a single cell type: In malignant epithelial cells: • Oxygen increases cycle entry towards 0.0007 1/min, with a Hill response, with half-max 21.5 mmHg and Hill power 4. • Pressure decreases cycle entry towards 0.0 1/min, with a linear response, with minimum threshold 0 and maximum threshold 0.5.• Dead increases debris secretion towards 1 1/min, with a Hill response, with half-max 0.1 and Hill power 4. Rule applies to dead cells.

Mathematical representation
With clearly defined behaviors and signals, and grammar to connect them, we can now uniquely map human-interpretable cell hypothesis statements onto mathematical expressions that make the grammar both human interpretable and computable.Moreover, our mathematical formulation allows new hypotheses to be directly added to models without modifying prior hypotheses, making our framing extensible and scalable as new knowledge is acquired.

Response functions
We use response functions to mathematically represent how a behavior varies with a signal towards its maximal (saturated) response.Our initial set of supported response functions are drawn from recurring forms in mathematical biology.In general, a response function  satisfies these properties: 1. (0) = 0.
In the absence of a signal, there is no response.2.  > () ≥ 0 Responses increase monotonically as the signals increases.

Linear response function
Many mathematical models use linear constitutive relations, so we define: Here,  9 and  " are (nonnegative) parameters of the response function governing where the response reaches its minimum and maximum values.For simplicity, we will call  9 and  " the minimum and maximum response thresholds, respectively.

Hill (sigmoidal) response function
Hill functions (or sigmoidal functions) are commonly used in pharmacodynamics and in computational models of cell responses to chemical signals.We define: and () = 0 if  < 0.
Here,  K25L and ℎ are (nonnegative) parameters governing where the response reaches half of its maximum value ( K25L ) and the steepness of the response (the Hill power ℎ).

Multivariate Hill (sigmoidal) response function
Sometimes, we may need to consider the impact of multiple signals  (with half-maximum parameters  K25L and Hill powers ) on a behavior.We therefore introduce a multivariate version: with the additional stipulation that we replace  ! by ( ! ) > = max( !, 0) as needed.It satisfies: Thus, the signals can contribute to a response individually and in combination, and when only one signal is supplied, the response behaves as the more common Hill response function.This functional form has the benefit that individual Hill response parameters do not need to be recalibrated as new arguments (i.e., new biological knowledge and rules) are added.See Fig. S1 (left).

Multivariate linear response function
Similarly, we may need to consider the impact of multiple signals  (with thresholds  9 and  " ) on a behavior.Motivated by the properties of the multivariate Hill function, we introduce a multivariate linear response function: with the additional stipulation that we replace  ! by ( ! ) > = max( !, 0) as needed.It satisfies: Thus, the signals can contribute to a response individually and in combination, and when only one signal is supplied, the response behaves as the simpler linear response function.This functional form has the benefit that individual linear response parameters do not need to be recalibrated as new arguments (i.e., new biological knowledge and rules) are added.See Fig. S1 (right).

Approximating linear responses with Hill response functions
When curating and aggregating prior knowledge, we may encounter linear relationships and seek to approximate them Hill response functions for a more unified framework.If a linear response function has thresholds  9 and _1, then we stipulate that: • The Hill response function reaches 90% of its peak value at  ) This can be satisfied when the Hill power is given by where 0 <  < 1 is a tolerance.(We set  = 0.1.)An Example is shown Fig. S2.We note that computationally, integer-valued Hill powers are at least tenfold faster in execution than non-integer powers, so it can be advantageous to round to the nearest whole number.

Approximating Hill responses with linear response functions
Similarly, when curating and aggregating prior knowledge, we may encounter Hill relationships and seek to approximate them as linear responses in a broader, unified framework.Moreover, this offers an opportunity for computational acceleration.If a Hill response function has half-max  K25L and Hill power ℎ, then we stipulate that: • The right threshold  ) is reached where the Hill response function reaches 90% of its peak value, and • The half-max and thresholds are related by  %&'( = This can be satisfied by setting: As before, 0 <  < 1 is a tolerance, which we set to  = 0.1.An Example is show in S2.

Single hypothesis statement: S increases B or S decreases B
For any biological hypothesis statement of the form "S increases B" or "S decreases B", the behavior B can be associated with a phenotypic behavioral parameter .The statement then takes the form of a linear interpolation in the nonlinear response function : () =  9 + ( N − p 9 )() = s1 − ()t ⋅  9 + () ⋅  J Here,  9 is the base level of the parameter (from the base behavioral phenotype),  J is the maximal (saturated) response, and () is a user-selected response function.

Example: oxygen-driven cycling
The statement "oxygen increases cycle entry" (with a "base" value of  9 = 0.001 hr O" and a maximum rate of  N = 0.042 hr O" ) takes the form If we use a linear response function as in prior work 6,7,10,13,21 with minimum and maximum thresholds at 5 mmHg and 38 mmHg, respectively, then the response function takes the form
If we use a Hill response function (with a half-max of 0.25, and Hill power of 3), the mathematical rule takes the form  9" = 0.001 + (0.0 − 0.001)

Multiple hypothesis statements: S1 increases B, S2 increases B
For any biological hypothesis statement of the form "S decreases B", the behavior B can be associated with a phenotypic behavioral parameter .The statement then takes the form: Here,  9 is the base level of the parameter (its value in the absence of any signals),  J is the maximal (saturated) response, and () is a user-selected response function.We use the multivariate Hill response function  J for ().

Example: oxygen-and hormone-dependent cycling
The statements "oxygen increases cycle entry" and "estrogen increases cycle entry" (with a "base" value of  9 = 0.001 hr O" and a maximum rate of  N = 0.042 hr O" ) takes the form

Multiple competing hypothesis statements: S1 increases B, S2 decreases B
For competing biological hypothesis statements of the form "S1 increases B" and "S2 decreases B", if the behavior has the associated phenotypic parameter , then we use the form: Here,  9 is the base level of the parameter (from the base behavioral phenotype in the absence of signals),  J is the maximal response to the promoting signal  " (with half-max  " * and Hill power ), and  : is the maximal response to the inhibiting signal  & (with half-max  & * and Hill power ).
We can write this more simply with: giving the overall response as a bilinear interpolation: In this formulation, the promoting signal  " pushes the behavior towards a "target" value that is then subject to inhibition by  & .Notice that: • If  * = 0, then  = 0 and ( ) ,  * ) behaves as the regular Hill response function ( ) ) for the promoting signal.

Example: "contradictory" observations on cell migration
Suppose a cell is migrating by chemotaxis towards a (scaled to be nondimensional) chemokine , where we observe: • Migration becomes less random and more directed as the chemokine concentration increases (it is more difficult to sense a chemical gradient in low concentrations) • Migration becomes more random and less directed for extremely high chemokine concentrations (high concentrations can saturate the cell's chemical receptors, making it difficult to sense the chemical gradient) In the cell behavior grammar, these statements become: •  increases migration bias o We will use a base migration bias of 0.01 and maximum bias of 1.
o If the effect reaches its maximum at  = 0.2, we set the linear response minimum and maximum thresholds to  + = 0 and  ) = 0.2.For the Hill response approximation, we set  %&'( = 0.1 and ℎ = 3 •  decreases migration bias o We will use a minimum migration bias of 0. o If the effect is first notable at  = 0.8 and increases towards its maximum at  = 1, we use minimum and maximum response thresholds of  + = 0.8 and  ) = 1 for a linear response.

General form: multiple promoting and inhibiting rules
Suppose we have the statements: Here, let  J be the maximum value of the parameter  (under the combined influence of the up-regulating signals  ), let  9 be its base value in the absence of signals, and let  : be its minimum value (under the combined influence of the down-regulating signals  ).
We define the total up response as: and the total down response as: Using these, we combine the overall response of the behavioral parameter as bilinear interpolation int eh nonlinear up-and down-responses  and : Notice that: • In the presence of up-regulating signals only, this reduces to the multivariate Hill response function  # (;  * , ).
• In the presence of down-regulating signals only, this reduces to the multivariate Hill response function  # (;  * , ) • Generally, the combined up-regulating signals sets a "target" value of the parameter, which can then be inhibited by the combined down-regulating signals.Note also that adding and removing individual rules to the form does not require alteration to the remaining rules.In this release, we use multivariate Hill response functions for clarity, but mixed linear and Hill responses could be used in the future.

Example:
We combine three hypothesis statements from prior examples that modulate cycle entry:

Reference PhysiCell implementation
We implemented support for the grammar in PhysiCell version 1.12.0 61 .Rules are imported at the start of a simulation using a compact CSV format (see supplementary information), parsed, and stored for each cell type.Rules are stored in a ruleset data structure, with a separate ruleset for each cell definition: 1.For each modulated behavior (with corresponding parameter ), store a set of rules: a.For each rule, we store: i.The "direction" of the response (increases/promotes or decreases/inhibits the behavior) ii.The signal used in the rule iii.The half-max, Hill parameter, and maximal value of the behavior of the parameter b.Based upon all the up-regulating rules, store the largest maximum parameter value  # c.Based upon all the down-regulating rules, store the lowest minimum parameter value  , To execute a set of rules (for a single cell): b.Query the base parameter value  + , the maximum parameter value  # , and the minimum parameter value  , .c. Compute the modulated parameter value via  = (1 − ) ⋅ [(1 − ) ⋅  + +  ⋅  # ] +  ⋅  , d. Update the phenotypic parameter  for the cell.At the start of each simulation step, each cell evaluates its individual rules (based on its current cell type) as noted above to set its current phenotype parameters, and then runs its standard phenotype processes.
We envision that other agent-based simulation frameworks (e.g., Chaste, CompuCell3D, Biocellion, Morpheus, and Tissue Simulation Toolkit) can independently implement this framework so long as they (1) implement the reference cell process models, (2) define cell types, (3) can generate compatible dictionaries of signals and behaviors, (4) can parse the rules, and (5) use these to modulate the cell behavioral parameters using the reference multivariate response functions as noted above.For an up-front development cost, simulation frameworks could advance reproducibility and support cross-model validation.

CELL
BEHAVIOR GRAMMAR: FULL DESCRIPTION Now that we have introduced the key signals and behaviors, we now describe our grammar to modulate cell behaviors in response to signals.Cell behavioral hypotheses are written as human-readable statements of the form: In [cell type T]: [signal S] [increases or decreases] [behavior B] [optional parameters 1 & 2].[optional statements].
When [base] and[effect] are not specified, parsers should assume a tenfold (one order-of-magnitude) change in behavior from the base parameter value: J =10  9 for increasing responses  : = 0.1  9 for decreasing responses Optional parameters 2: response form and parameters Any rule can specify the form of its response function (currently linear or Hill) by appending: with a [linear or Hill] response.If the rule further specifies the parameters of the response, it takes the form: with a [linear or Hill] response, with [parameters].For linear responses, [parameters] takes the form: with minimum threshold [value] and maximum threshold [value].For Hill responses, [parameters] takes the form: with half-max [value] and Hill power [value].

Fig. S1 :
Fig. S1: Multivariate response functions.Left: A multivariate Hill response function, where  ! has half-max 2 and Hill power 8, and  " has half-max 1 and Hill power 2. Half-maxes are plotted on each variable as blue dashed lines, while the half-max across the multivariate function is plotted as a red contour.( ! ) and ( " ) are plotted as black curves along the respective axes.Right: A multivariate linear response function, where  ! has min and max thresholds 0.75 and 3.75, and  " has thresholds 0 and 2. Half-maxes are plotted on each variable as blue dashed lines, while the half-max across the multivariate function is plotted as a red contour.( ! ) and ( " ) are plotted as black curves along the respective axes.

𝑟 9
If we use a Hill response function (with a half-max of 21.5 mmHg, and Hill power of 4), the mathematical rule takes the form  9" = 0.001 + (0.042 − 0.001) (pO & ) P 21.5 P + (pO & ) P .Fig. S2: Converting between linear and Hill response functions.Left: Approximating a linear response function (black dotted curve) with a Hill response function (red).Right: Approximating a Hill response function (black dotted) with a linear response function (red).

Fig. S3 :
Fig. S3: Sample response functions.Left: The behavioral response function for the statements "oxygen increases cycle entry" and "pressure decreases cycle entry", using the multivariate Hill response function.Right: A non-monotonic response function from the statements "c increases migration bias" (for lower values of c) and "c decreases migration bias" (for higher values of c), showing both linear and Hill response constructions.

1 .
For each modulated behavior (with corresponding parameter  ): a.For each rule: i. Sample the corresponding signal  ii.Query its Hill parameter ℎ and half-max  %&'( iii.Compute S / / #$%& T 0 iv.If the rule increases (promotes) the behavior, add the contribution from iii to  in the generalized multivariate response function.Otherwise, add the contribution to .

Example 1 :
Model the progression of hypoxia in a metastatic tumor Example 2: Forecast locations of tumor cell invasion from an initial state defining tumor cell phenotypes with spatial transcriptomics Example 4: T cell activation, expansion, and exhaustion in a diverse tumor microenvironment Example 5: Using experimental insight to model combination immune-targeted therapies in the pancreatic cancer microenvironment

TABLE 2 .
Death terms in the behavior dictionary, along with corresponding controllable parameters in the PhysiCell reference implementation.

Table 8 .
A list of symbols in our vocabulary of signals that can be used to drive behavioral or state changes in cells.Here, X is any diffusible substrate, Y is a cell type, and Z is a customized cell parameter (or variable).