Skip to main content
  • Research article
  • Open access
  • Published:

Generation of multicellular spatiotemporal models of population dynamics from ordinary differential equations, with applications in viral infection

Abstract

Background

The biophysics of an organism span multiple scales from subcellular to organismal and include processes characterized by spatial properties, such as the diffusion of molecules, cell migration, and flow of intravenous fluids. Mathematical biology seeks to explain biophysical processes in mathematical terms at, and across, all relevant spatial and temporal scales, through the generation of representative models. While non-spatial, ordinary differential equation (ODE) models are often used and readily calibrated to experimental data, they do not explicitly represent the spatial and stochastic features of a biological system, limiting their insights and applications. However, spatial models describing biological systems with spatial information are mathematically complex and computationally expensive, which limits the ability to calibrate and deploy them and highlights the need for simpler methods able to model the spatial features of biological systems.

Results

In this work, we develop a formal method for deriving cell-based, spatial, multicellular models from ODE models of population dynamics in biological systems, and vice versa. We provide examples of generating spatiotemporal, multicellular models from ODE models of viral infection and immune response. In these models, the determinants of agreement of spatial and non-spatial models are the degree of spatial heterogeneity in viral production and rates of extracellular viral diffusion and decay. We show how ODE model parameters can implicitly represent spatial parameters, and cell-based spatial models can generate uncertain predictions through sensitivity to stochastic cellular events, which is not a feature of ODE models. Using our method, we can test ODE models in a multicellular, spatial context and translate information to and from non-spatial and spatial models, which help to employ spatiotemporal multicellular models using calibrated ODE model parameters. We additionally investigate objects and processes implicitly represented by ODE model terms and parameters and improve the reproducibility of spatial, stochastic models.

Conclusion

We developed and demonstrate a method for generating spatiotemporal, multicellular models from non-spatial population dynamics models of multicellular systems. We envision employing our method to generate new ODE model terms from spatiotemporal and multicellular models, recast popular ODE models on a cellular basis, and generate better models for critical applications where spatial and stochastic features affect outcomes.

Background

The function and form of biological systems involve biophysical mechanisms across all scales, from the subcellular to the organismal [1]. Processes and properties observed at a particular scale emerge from, and affect, complex processes at a finer scale, such as the emergent polarization and translocation of a cell from complex subcellular reaction kinetics when exposed to a chemotactic stimulus [2]. Such hierarchical organization can be established to relate the roles of individual biochemical components in signaling and cytoskeletal networks at the subcellular level to the emergence of tissue-level force generation, cell migration and tissue function, and shape changes [3]. Between these two scales lies the cell, which has been argued to provide a natural level of abstraction for modeling development [4]. The role played by individual cells and their interactions is apparent in development and disease, including tip cells in collective cell migration [5], diversity in production of virions during viral infection [6], and cancer initiation by stem cells [7].

Each cell occupies a volume in space, the biochemical processes in which are well known to experience effects of spatial mechanisms like volume exclusion of macromolecules, the cytoskeleton, and organelles, resulting in macromolecular crowding that causes highly non-mass-action reaction rates [8] and anomalous diffusion [9]. Recent computational models have sought to describe the spatiotemporal effects on the kinetics of elementary reactions [10], as well as to relate discrete and continuous mathematical descriptions of spatially resolved subcellular reaction kinetics [11]. Likewise, models of spatiotemporal multicellular systems in development and disease have shown the significance of spatial mechanisms like diffusive transport and the shape and position of individual cells in angiogenesis [12], polycystic kidney disease [13], and spheroid fusion [14]. Simulations of viral infection have demonstrated the non-negligible effects of the well-mixed assumptions commonly employed when modeling viral infection and immune response using population dynamics, like the neglect of the initial distribution of infected cells in susceptible tissue [15].

The ability to derive cell-based, spatiotemporal models from ordinary differential equation (ODE) models would enhance the utility of both types of models. Cell-based, spatiotemporal models can explicitly describe cellular and spatial mechanisms neglected by ODE models that affect the emergent dynamics and properties of multicellular systems, such as the influence of dynamic aggregate shape on diffusion-limited growth dynamics [16] and individual infected cells on the progression of viral infection [17]. Likewise, ODE models can inform cell-based, spatiotemporal models with efficient parameter fitting to experimental data, and can appropriately describe dynamics at coarser scales and distant locales with respect to a particular multicellular domain of interest (e.g., the population dynamics of a lymph node when explicitly modeling local viral infection). One such example is the approach of Murray and Goyal to derive discrete stochastic dynamics from continuous dynamical descriptions using the Poisson distribution in their multiscale modeling work on hepatitis B virus infection [18]. Likewise Figueredo et al. compared derived representations for mechanisms associated with early-stage cancer using non-spatial agent-based, ODE and stochastic differential equation modeling approaches and demonstrated the feasibility of generating equivalent mechanistic models [19]. However, to our knowledge, no well-defined general formalism describes systematic translation of models to the cellular scale from coarser scales at which spatially homogeneous, population dynamics models using ODEs appropriately describe a biological system. In the very least, the lack of consistent translation of model terms and parameters between spatial and non-spatial models severely inhibits the potential to apply the vast amount of available information and resources in non-spatial modeling, such as those available in BioModels [20], to spatial contexts, and to share validated, reproducible spatial models [21].

In this work, we develop a method for generating spatial, cell-based models from homogeneous, ODE-based models of biological systems to address critically important questions such as

  1. 1.

    Is the behavior of an ensemble of N cells with mean parameter x in mechanism y the same as N times the behavior of a single cell with the same mechanism and mean parameter?

  2. 2.

    Is the behavior of a biological system in which a species realistically diffuses qualitatively different from one in which the species diffuses infinitely fast?

As such, we develop our method under the premise that a spatial, cell-based model applied to an ensemble of cells reproduces the ODE model from which it was generated in the limit of well-mixed conditions. We describe the general formalism of our method with regard to continuous modeling of cell populations and soluble signals, and spatiotemporal, multicellular, and cell-based mechanisms (e.g., diffusion, contact-mediated interactions, discrete cell type transition probabilities). We apply our method to generate spatial, cell-based models from two non-spatial models of viral infection and host-pathogen interaction based on our recent work in viral infection and immune response modeling [17] to show critically important aspects of spatially homogeneous vs. heterogeneous models, and to demonstrate how the novel ability through our method to generate explicit models of implicit mechanisms in a homogeneous model can be used to generate new hypotheses about a biological system. We call our method cellularization, which is to mean that cellularizing a homogeneous model generates an explicit representation of cellular-level mechanisms implicitly represented by the homogeneous model.

Results

Spread of viral pathogens during infection is well known to be an inherently spatial process, though the spatiotemporal details of viral spread have not been well studied [22]. Recent computational modeling and simulation have shown the critical importance of considering spatial spread of infection when studying immune response mechanisms like autocrine and paracrine interferon signaling [23], while recent single-cell transcriptomics has shown the significance of cellular heterogeneity in the innate immune response to influenza A virus [24]. In this section, we demonstrate the cellularization of ODE models of viral infection and immune response according to the formalism defined in “Conclusion.”

In this work, we developed a method for generating spatial, multicellular models of biological systems from non-spatial models, and vice versa, which we call cellularization. We demonstrate using our method by cellularizing non-spatial models of viral infection and host-pathogen interaction. Using these cellularized models, we quantitatively showed how spatial mechanisms implicitly represented in non-spatial models can exhibit significant effects on emergent dynamics when explicitly modeled. Variations in related non-spatial model parameters emerged from moderate cases of varying spatial mechanisms like rate of diffusive mass transport of virus, while extreme cases generated emergent dynamics inconsistently with those described by the non-spatial model. We describe the responsible mechanisms for extreme disagreement between homogeneous and cellular models, specifically concerning the limitations of describing discrete biological objects and processes using continuous descriptions.

Methods and deployment of their analogous spatial models in simulation. We present ODE and spatial model simulation results while investigating emergent effects related to spatial, cell-based, and stochastic aspects introduced by the spatial models. Both cellularized ODE models were generated ad hoc for the purposes of this work according to published models of viral infection and immune response [25, 26], the parameters of which were selected only for conveniences of demonstrating their spatial analogues. As such, the simulation results in this work per se should not be regarded as relevant to modeling a particular virus or patient scenario. Rather, simulation results should be surveyed within the scope of this work while considering the possible applications of cellularization to existing models of a particular virus or patient scenario.

Two-dimensional infectivity

The first cellularized ODE model in this work consists of the interactions of a population of cells and an extracellular virus using mass action. The virus infects uninfected cells U with an infection rate β, which causes uninfected cells to join an infected cell population I1. Infected cells become virus-releasing cells I2 at a rate k, at which point they release unitless virus V at a rate p, modeling an eclipse phase. Virus-releasing cells die and join a dead cell population A at a rate d, modeling virally induced cell death. The virus decays at a rate c,

$$ \left\{\begin{array}{c}\frac{dU}{dt}=-\beta UV\\ {}\frac{d{I}_1}{dt}=\beta UV-k{I}_1\\ {}\frac{d{I}_2}{dt}=k{I}_1-d{I}_2\\ {}\frac{dA}{dt}=d{I}_2\\ {}\frac{dV}{dt}=p{I}_2- cV\end{array}\right. $$
(1)

To build a spatial model of local infection using the configuration described in “Implementation details,” Eq. 1 can be rescaled to a simulated tissue patch and considered as global measurements of the spatial domain. Knowing the total number of epithelial cells in both the ODE and spatial models, the global scaling factor η can be calculated from Eq. 12 and applied to generate a form that describes the quantities in Eq. 1 but at the size of the spatial domain,

$$ \left\{\begin{array}{c}\frac{du}{dt}=-\frac{\beta }{\eta } uv\\ {}\frac{d{i}_1}{dt}=\frac{\beta }{\eta } uv-k{i}_1\\ {}\frac{d{i}_2}{dt}=k{i}_1-d{i}_2\\ {}\frac{da}{dt}=d{i}_2\\ {}\frac{dv}{dt}=p{i}_2- cv\end{array}\right.. $$
(2)

as in “Conclusion.”

In this work, we developed a method for generating spatial, multicellular models of biological systems from non-spatial models, and vice versa, which we call cellularization. We demonstrate using our method by cellularizing non-spatial models of viral infection and host-pathogen interaction. Using these cellularized models, we quantitatively showed how spatial mechanisms implicitly represented in non-spatial models can exhibit significant effects on emergent dynamics when explicitly modeled. Variations in related non-spatial model parameters emerged from moderate cases of varying spatial mechanisms like rate of diffusive mass transport of virus, while extreme cases generated emergent dynamics inconsistently with those described by the non-spatial model. We describe the responsible mechanisms for extreme disagreement between homogeneous and cellular models, specifically concerning the limitations of describing discrete biological objects and processes using continuous descriptions.

Methods, here u measures the same quantity as U, i1 measures the same quantity as I1, and so on, but at the scale of the spatial domain. The ODE model describes the processes of virus-mediated infection an uninfected cell (\( \hat{U}\to \hat{I_1} \)), an infected cell becoming a virus-releasing cell (\( \hat{I_1}\to \hat{I_2} \)), and a virus-releasing cell dying due to the effects of the virus (\( \hat{I_2}\to \hat{A} \)). Spatial analogues of these processes can be formulated using Eq. 29 for a diffusive extracellular virus field \( \overset{\sim }{v} \) of the form from Eq. 22,

$$ \left\{\begin{array}{c}\Pr \left(\tau \left(s,t+\Delta t\right)=\hat{I_1}|\tau \left(s,t\right)=\hat{U}\right)=1-{e}^{-\frac{\beta }{\theta}\overline{v}(s)\Delta t}\\ {}\Pr \left(\tau \left(s,t+\Delta t\right)=\hat{I_2}|\tau \left(s,t\right)=\hat{I_1}\right)=1-{e}^{-k\Delta t}\\ {}\Pr \left(\tau \left(s,t+\Delta t\right)=\hat{A}|\tau \left(s,t\right)=\hat{I_2}\right)=1-{e}^{-d\Delta t}\\ {}{\partial}_t\overset{\sim }{v}={D}_V{\partial}_i^2\overset{\sim }{v}-c\overset{\sim }{v}+\overset{\sim }{w}\\ {}\overset{\sim }{w}\left(x,t\right)=\left\{\begin{array}{cc}\frac{p}{\left\Vert \mathcal{V}\left(\sigma \left(x,t\right),t\right)\right\Vert }& \tau \left(\sigma \left(x,t\right),t\right)=\hat{I_2}\\ {}0& \mathrm{otherwise}\end{array}\right.\end{array}\right.. $$
(3)

Cellular measurements of extracellular virus \( \overline{v} \) are calculated from \( \overset{\sim }{v} \) according to Eq. 19. Unless otherwise specified, all simulations were performed with the parameter values in Tables 1 and 3. Simulations were initialized with fractions of initially infected cells in the range 1% [28] to 10% [29], the distributions of which were random in the spatial model, and with a single initially infected cell, which was placed at the center of the epithelial sheet in the spatial model [17].

Table 1 Model parameters used in simulations of infection unless otherwise specified

The distribution of infected and virus-releasing cells demonstrated notable spatial features early in simulation, and especially for fewer initially infected cells (Fig. 1). Groups of infected and virus-releasing epithelial cells typically formed near sites of initial infection, which generated prominent gradients in the distribution of extracellular virus. However, some sites of significant infection were also observed early where cells were not initially infected (e.g., Fig. 1, 0.01 initial infection fraction, 3 days in simulation time), which produced notable gradients in the extracellular virus field. These sites of infection coalesced into infection throughout the simulated tissue patch, the earliest infected cells of which often, but not always, corresponded to outgrowth of dead cells, and produced elevated levels of extracellular virus throughout the spatial domain. Without any antiviral strategy, nearly all cells died by 2 weeks of simulation time for all initial infection fractions.

Fig. 1
figure 1

Spatial model results of viral infection in a two-dimensional, epithelial sheet. Results shown for 1% (top), 5% (middle), and 10% (bottom) initially infected cells at 0, 3, 4, 5, 7, and 14 days in simulation time. Epithelial cells shown as blue when susceptible, green when infected, red when virus-releasing, and black when dead. Lower-right color bar shows levels in the virus field, from blue (0) to red (0.2)

To test the validity of cellularization and the effects of spatial mechanisms introduced to the ODE model, scalar measurements of spatial model results were made of all cell type populations and the extracellular virus and compared to results from the scaled ODE model (Eq. 2). By inspection, the spatial model and employed model parameters generated spatiotemporal dynamics consistently with the original ODE model (Fig. 2). Among the 100 simulation replicas of the spatial model for each initial infection fraction, an initial infection fraction of 0.01 generated results with the most significant differences compared to those from the scaled ODE model (worst and best QRF 55% and 84%, infected cells and virus-releasing cells, respectively; worst and best NMAE 30% and 4%, dead cells and susceptible cells, respectively). In the case of an initial infection fraction of 0.01, some simulation replicas produced infection dynamics with a delay, where all measures of progression of infection occurred at the same magnitude and with the same dynamical features, but slightly later than predicted by the ODE model, creating a rightward skew in replica results. Differences between models decreased with increase initial infection fraction, where 0.1 initial infection fraction produced the best agreement (worst and three best QRF 70% and 100%, susceptible cells and the three variables dead cells, virus-releasing cells and virus, respectively; worst and best NMAE 20% and 1%, susceptible cells and dead cells, respectively), and 0.05 produced more, but still seemingly acceptable by inspection, differences (worst and best QRF 64% and 100%, susceptible cells and dead cells, respectively; worst and best NMAE 22% and 1%, susceptible cells and dead cells, respectively).

Fig. 2
figure 2

Scalar results from 100 replicas of the spatial model of viral infection. Results shown for 1% (left), 5% (center), and 10% (right) initially infected cells. Medians shown as black lines. 0th to 100th quantiles shaded as blue, 10th to 90th as orange, and 25th to 75th as light blue. ODE model results shown as red dashed lines

Spatial effects were even more prominent for simulation replicas with one initially infected cell (Fig. 3). Spread of infection exhibited outward spread from the initial site of infection, but with virus-releasing cells far from the initial site of infection even by 3 days of simulation time (Fig. 3A). Such spread of infection at seemingly random locations was especially noticeable by 7 days of simulation, where the distributions of dead cells were comparably distributed as those for an initial infection fraction of 0.1, though the degree of death occurred later in simulation. The observed delay due to lower initial infection was increasingly notable for simulation replicas with one initially infected cell, where peak viral load occurred for some simulation replicas on the order of days later than when predicted by the ODE model (Fig. 3B). For some simulation replicas with one initially infected cell, progression of infection did not occur at all. Replicas for one initially infected cell produced significantly more variability in results, which produced some better QRF measurements compared to initial infection fractions from 0.01 to 0.1 (worst and best QRF 80% and 92%, susceptible cells and virus, respectively), but significantly worse NMAE measurements (worst and best NMAE 68% and 20%, infected cells and dead cells, respectively).

Fig. 3
figure 3

Spatial and ODE model results of viral infection in a two-dimensional, epithelial sheet with one initially infected cell. A Spatial model results of one simulation replica at 0, 3, 4, 5, 7, and 14 days in simulation time. Epithelial cells and extracellular virus shown as in Fig. 1. B Scalar results from 100 replicas of the spatial model. Medians shown as black lines. 0th to 100th quantiles shaded as blue, 10th to 90th as orange, and 25th to 75th as light blue. ODE model results shown as red dashed lines

Effects of diffusivity

Cellularization of ODE models like Eq. 2 presents the ability to interrogate the effects of mass transport that are implicitly represented in a homogeneous model. To test the effects of extracellular diffusion of infection dynamics, the virus diffusion coefficient was swept for all aforementioned initial infection conditions with logarithmic variations over a range of a reduction by a factor of 1000, to an increase by a factor of 10, from the value in Table 1 (Fig. 4). Decreases in virus diffusivity generally led to reductions in the rate of infection of the tissue patch and greater departure from the dynamics described by the ODE model. As in results shown in Fig. 3, no infection occurred for some replicas with one initially infected cell and all diffusion coefficients less than 1.0 μm2/s. In cases with the minimum considered virus diffusion coefficient (i.e., 0.0001 μm2/s), the number of susceptible cells tended towards a non-zero final value, and variations with very small virus diffusion coefficients and smaller initial infection fractions generated replicas in strong disagreement with the ODE model (e.g., for 0.0005 μm2/s and one initially infected cell, worst QRF 2%, susceptible cells). For initial infection fractions of 0.05 and 0.1, the spatial and ODE models agreed for virus diffusion coefficients as low as 0.01 μm2/s, whereas the same was true for an initial infection fraction of 0.01 and virus diffusion coefficient as low as 0.05 μm2/s. With increasing initial infection fraction and virus diffusion coefficient, replicas tended towards the ODE model (e.g., for virus diffusion coefficient of 1 μm2/s and all initial infection fractions of at least 0.01, QRF was 1 for all variables and worst NMAE was 7% for susceptible cells and 0.01 initial infection fraction).

Fig. 4
figure 4

Infection model results while varying spatial model parameters. A Number of susceptible cells from 100 replicas of the spatial model of viral infection while varying initial infection fraction and virus diffusion coefficient. Medians shown as black lines. 0th to 100th quantiles shaded as blue, 10th to 90th as orange, and 25th to 75th as light blue. ODE model results shown as red dashed lines. ODE model results with best-fit infectivity shown as magenta dotted lines. B Efficiency of infectivity (top), measured as the ratio of fitted to original ODE model infectivity, and score of fit (bottom), measured as the mean of the NMAE for all variables of the fitted ODE model, for all spatial model parameter variations. Results for efficiency of infectivity are shaded from lowest value (red) to highest value (green), and for score of fit from zero (green) to one (red)

We tested whether decreases in virus diffusion coefficient in the spatial model correspond to decreases in infection rate (i.e., the parameter β) by fitting the ODE model to spatial model results while varying the infection rate. Fitting was performed using results from all simulation replicas per parameter set and initial conditions when infection occurred. The error during fitting was calculated as the squared difference between the fitted ODE model results and spatial model results for all cell type populations and the viral load at all time points. We calculated an efficiency of infectivity as the ratio of the infectivity as fitted to the spatial model results to the infectivity of the original ODE model. We also calculated a score of fit to characterize the overall similarity of the fitted ODE model and spatial model results as the mean NMAE of all variables of the fitted ODE model. Reduced infection rates were found that reproduced spatial model results using the ODE model for all cell type populations and viral load (Fig. 5) for many virus diffusion coefficients, with a range of efficiency of infectivity between 6.6% (single initially infected cell, virus diffusion coefficient of 0.0001 μm2/s, score of fit of 0.262) and 99.2% (initial infection fraction of 0.05, virus diffusion coefficient of 1.0 μm2/s, score of fit of 0.026). The ODE model better reproduced spatial model results with reduced infection rates for greater virus diffusion coefficients and initial infection fractions (e.g., for fitted ODE model to 0.01 μm2/s virus diffusion coefficient and 0.05 initial infection fraction, worst and three best QRF: 90% and 100%, infected cells and the three variables virus-releasing cells, dead cells and virus, respectively; worst and best NMAE: 25% and 0.9%, susceptible cells and dead cells, respectively). Some virus diffusion coefficients generated dynamics that were not consistent with the ODE model and reduced infection rates, such as virus diffusion coefficients less than 0.005 μm2/s for an initial infection fraction of 0.01, and less than 0.001 μm2/s for initial infection fractions of 0.05 and 0.1 (e.g., for fitted ODE model to 0.0001 μm2/s virus diffusion coefficient and 0.1 initial infection fraction, worst QRF of 27% for virus-releasing cells). Otherwise, the ODE model could be well-fitted to spatial model results with decreasing virus diffusion coefficient by reducing the ODE model infection rate (e.g., all QRF measurements of fitted ODE models were at least 94% for all virus diffusion coefficients of at least 0.1 μm2/s).

Fig. 5
figure 5

Viral load from 100 replicas of the spatial model of viral infection while varying initial infection fraction and virus diffusion coefficient. Medians shown as black lines. 0th to 100th quantiles shaded as blue, 10th to 90th as orange, and 25th to 75th as light blue. ODE model results shown as red dashed lines. ODE model results with best-fit infectivity shown as magenta dotted lines

Immune response modeling

The second cellularized ODE model in this work adds a compartmental model of immune cell proliferation and recruitment to the ODE model of viral infection in “Two-dimensional infectivity,” which presents a number of spatially resolved mechanisms well suited for studying using cellularization. In this model, we neglect most of the intricacies of autocrine and paracrine signaling involved in the innate and adaptive immune responses that, though critical to an effective immune response in vivo and candidates for applications of cellularization in model-based research, need not be explicitly modeled for the purposes of this demonstration (for an example of detailed modeling of the immune response using ODEs, see [25]). Infected and virus-releasing cells release a local cytokine C at a rate pC, which decays at a rate cC. The local cytokine transports to cytokine in a lymph node compartment CL at a rate kC, which decays at a rate ccl. The lymph node cytokine induces proliferation of an immune cell population in the lymph node compartment EL at a rate pel, which also undergoes saturated proliferation and transports to a local immune cell population E at a rate ke. Each local immune cell kills virus-releasing cells at a rate dei2 [25], and the local immune cell population decays at a rate dE,

$$ \left\{\begin{array}{c}\frac{dU}{dt}=-\beta UV\\ {}\frac{d{I}_1}{dt}=\beta UV-k{I}_1\\ {}\frac{d{I}_2}{dt}=k{I}_1-\left(d+{d}_{ei2}E\right){I}_2\\ {}\frac{dA}{dt}=\left(d+{d}_{ei2}E\right){I}_2\\ {}\frac{dV}{dt}=p{I}_2- cV\\ {}\frac{dC}{dt}={p}_C\left({I}_1+{I}_2\right)-{c}_CC-{k}_CC\\ {}\frac{d{C}_L}{dt}={k}_CC-{c}_{cl}{C}_L\\ {}\frac{dE}{dt}={k}_e{E}_L-{d}_EE\\ {}\frac{d{E}_L}{dt}={p}_{el}{C}_L+{r}_{el}{C}_L{E}_L\frac{K_{el}}{K_{el}+{E}_L}-{k}_e{E}_L\end{array}\right. $$
(4)

Continuing the construction of a spatial model of a local patch of infection as described in “Two-dimensional infectivity,” Eq. 4 can be scaled in the same manner as performed on Eq. 1 using the total number of epithelial cells in the spatial and ODE models,

$$ \left\{\begin{array}{c}\frac{du}{dt}=-\frac{\beta }{\eta } uv\\ {}\frac{d{i}_1}{dt}=\frac{\beta }{\eta } uv-k{i}_1\\ {}\frac{d{i}_2}{dt}=k{i}_1-\left(d+\frac{d_{ei2}}{\eta }e\right){i}_2\\ {}\frac{da}{dt}=\left(d+\frac{d_{ei2}}{\eta }e\right){i}_2\\ {}\frac{dv}{dt}=p{i}_2- cv\\ {}\frac{dc}{dt}={p}_C\left({i}_1+{i}_2\right)-{c}_Cc-{k}_Cc\\ {}\frac{d{c}_L}{dt}={k}_Cc-{c}_{cl}{c}_L\\ {}\frac{de}{dt}={k}_e{e}_L-{d}_Ee\\ {}\frac{d{e}_L}{dt}={p}_{el}{c}_L+{r}_{el}{c}_L{e}_L\frac{K_{el}}{\eta {K}_{el}+{e}_L}-{k}_e{e}_L\end{array}\right. $$
(5)

Note that rescaling the cytokine and immune cell population in the lymph node compartment is optional and could instead remain at the scale of the original ODE model (e.g., eL = ηEL throughout). In the spatial model, local cytokine and each local immune cell are spatially modeled, whereas the lymph node cytokine and immune cell population remain as ODEs. Contact-mediated killing, which involves a number of spatially resolved properties (e.g., position) and dynamics (e.g., motility) that are implicitly represented in the ODE model, is cellularized for each virus-releasing cell using Eq. 34 for \( \hat{E}- \) type (i.e., local immune cell type) cell total contact area AE, total number of epithelial cells Ne, and total contact area between a cell s and \( \hat{E} \) cells Aσ, s. As such, Eq. 4 can then be written according to the aforementioned multiscale structure,

$$ \left\{\begin{array}{c}\Pr \left(\tau \left(s,t+\Delta t\right)=\hat{I_1}|\tau \left(s,t\right)=\hat{U}\right)=1-{e}^{-\frac{\beta }{\theta}\overline{v}\left(s,t\right)\Delta t}\\ {}\Pr \left(\tau \left(s,t+\Delta t\right)=\hat{I_2}|\tau \left(s,t\right)=\hat{I_1}\right)=1-{e}^{-k\Delta t}\\ {}\Pr \left(\tau \left(s,t+\Delta t\right)=\hat{A}|\tau \left(s,t\right)=\hat{I_2}\right)=1-{e}^{-\left(d+\frac{d_{ei2}{N}_e}{\eta {A}_E}{A}_{s,E}\right)\Delta t}\\ {}{\partial}_t\overset{\sim }{v}={D}_V{\partial}_i^2\overset{\sim }{v}-c\overset{\sim }{v}+{\overset{\sim }{w}}_V\\ {}{\partial}_t\overset{\sim }{c}={D}_C{\partial}_i^2\overset{\sim }{c}-\left({c}_C+{k}_C\right)\overset{\sim }{c}+{\overset{\sim }{w}}_C\\ {}{\overset{\sim }{w}}_V\left(x,t\right)=\left\{\begin{array}{cc}\frac{p}{\left\Vert \mathcal{V}\left(\sigma \left(x,t\right),t\right)\right\Vert }& \tau \left(\sigma \left(x,t\right),t\right)=\hat{I_2}\\ {}0& \mathrm{otherwise}\end{array}\right.\\ {}{\overset{\sim }{w}}_C\left(x,t\right)=\left\{\begin{array}{cc}\frac{p_C}{\left\Vert \mathcal{V}\left(\sigma \left(x,t\right),t\right)\right\Vert }& \tau \left(\sigma \left(x,t\right),t\right)\in \left\{\hat{I_1},\hat{I_2}\right\}\\ {}0& \mathrm{otherwise}\end{array}\right.\\ {}\Pr \left(\mathrm{add}\ k\ \hat{E}-\mathrm{type}\ \mathrm{cells}\right)=1-{e}^{-{k}_e{e}_L\Delta t}\sum \limits_{0\le n\le k}\frac{{\left({k}_e{e}_L\Delta t\right)}^n}{n!}\\ {}\Pr \left(\mathrm{remove}\ s|\tau \left(s,t\right)=\hat{E}\right)=1-{e}^{-{d}_E\Delta t}\end{array}\right., $$
(6)

where the rate equations for eL and cL are the same as in Eq. 5. All other equations are described in “Two-dimensional infectivity.” Cellularization of an ODE model like Eq. 4 also necessitates an explicit description of where and how incoming and outgoing cells enter and exit a local spatial domain, the details of which are neglected in a homogeneous model. For the cellularization in this work, we designated the spatial domain boundary above the epithelial sheet was assigned as the inflow boundary for local immune cells. Outgoing cells were immediately removed from the spatial domain and replaced with medium. Unless otherwise specified, all simulations were performed with the parameter values in Tables 1, 2 and 3.

Table 2 Model parameters used in simulations of infection and immune response unless otherwise specified
Table 3 Parameter values used in all simulations

The immune response prevented complete spread of infection throughout the spatial domain in all 100 simulation replicas of 0.01, 0.05, and 0.1 initial infection fraction, and decreasingly so with increasing initial infection fraction (Fig. 6). After 2 days of simulation time, infected and virus-releasing cells were more prominently distributed throughout the spatial domain compared to local immune cells. However, a strong immune response recruited sufficient immune cells to nearly cover the entire epithelial sheet by days 7, 6, and 5 for initial infection fractions 0.01, 0.05, and 0.1, respectively, resulting in significant killing of virus-releasing cells and prevention of further infection. By 2 weeks of simulation time for all initial conditions and simulation replicas, virus and cytokine levels were near zero, most immune cells had left the spatial domain and the epithelial sheet consisted of susceptible cells with significant distributions of dead cells. Initial infection fractions of 0.01, 0.05, and 0.1 resulted in final fractions of dead cells of around 0.5, 0.625, and 0.75, respectively.

Fig. 6
figure 6

Spatial model results of viral infection and immune response in a quasi-two-dimensional, epithelial sheet. Results shown for 1% (top), 5% (middle), and 10% (bottom) initially infected cells at 0, 2, 3, 4, 7, and 14 days in simulation time. Epithelial cells shown as in Fig. 1. Immune cells shown as dark red. Lower-right color bar shows levels in the virus and cytokine fields, from blue (0) to red (0.05)

As with the model of viral infection described in “Two-dimensional infectivity,” the spatial model of viral infection and immune response generated spatiotemporal dynamics consistently with the original ODE model using the employed model parameters (Fig. 7). Simulation replicas of the spatial model produced nearly the same number of susceptible cells at the end of simulation as the ODE model for 0.01 initial infection fraction, while replicas of the spatial model produced slightly fewer final susceptible cells and more final dead cells than the ODE model for initial infection fractions of 0.05 and 0.1. For an initial infection fraction of 0.01, peak viral load, local cytokine, and lymph node cytokine in replicas using the spatial model occurred slightly later than for the ODE model (worst and best QRF 30% and 88%, lymph node cytokine and susceptible cells, respectively; worst and best NMAE 40% and 5%, virus-releasing cells and susceptible cells, respectively), while these metrics were approximately the same between the spatial and ODE models for initial infection fractions of 0.05 and 0.1 (for 0.05 initial infection fraction, worst and best QRF 25% and 70%, dead cells and infected cells, respectively; worst and best NMAE 30% and 10%, virus-releasing cells and dead cells, respectively; for 0.1 initial infection fraction, worst and best QRF 16% and 70%, dead cells and infected cells, respectively; worst and two best NMAE 32% and 10%, virus-releasing cells and the two variables dead cells and local immune cells, respectively). The number of infected, virus-releasing, and lymph immune cells significantly varied for replicas of the spatial model and an initial fraction of 0.01, while for initial infection fractions of 0.05 and 0.1 the maximum number of infected, virus-releasing, and lymph immune cells were slightly greater than those of the ODE model. However, spatial model replicas produced slightly fewer local immune cells than the ODE model before peak viral load and then tended towards the same values as the ODE model thereafter.

Fig. 7
figure 7

Scalar results from 100 replicas of the spatial model of viral infection and immune response. Results shown for 1% (left), 5% (center), and 10% (right) initially infected cells. Medians shown as black lines. 0th to 100th quantiles shaded as blue, 10th to 90th as orange, and 25th to 75th as light blue. ODE model results shown as red dashed lines

In some simulation replicas of the spatial model, viral load suddenly increased to appreciable values even days after viral load decreased below a value of 1, and then decreased again to negligible values. Inspection of simulation data showed that latent extracellular virus infected susceptible cells after the initial progression of infection in some replicas of the spatial model. We observed this event, which cannot occur using the ODE model, in previous work modeling viral infection and called it “recursion” [17].

To demonstrate the utility of cellularization to generate explicit spatial models of underlying biophysical mechanisms and use them to formulate biological hypotheses, we tested the effects of spatial mechanisms introduced during development of the spatial analogue to the ODE model of infection and immune response. Specifically, we performed a parameter sweep of the sampling fraction by which immune cells are seeded and the chemotaxis parameter that describes how sensitively they chemotax along gradients of local cytokine for an initial infection fraction of 0.1. We varied the sampling fraction from 2.5 × 10−5 (i.e., equivalent to seeding in random locations) to 1.0 (i.e., at the maximum value of local cytokine in all available sites), and the chemotaxis parameter from 104 (i.e., virtually no chemotaxis) to 106 (i.e., just below permitting annihilation of cells by excessive chemotactic forces) using a logarithmic scale. We measured the normalized root mean squared error (NRMSE) of results from each of 100 simulation replicas as forecasting results from the ODE model. NRMSE calculations were made for results in intervals of ten simulation steps (i.e., 50 min) for each parameter set, the summation of which we used to quantify how well the spatial model represented the constituent components of the ODE model using each parameter set.

As determined by the final number of susceptible cells in simulation replicas, we found that the most effective immune response strategy was one of random seeding and strong chemotaxis, and the least effective strategy was one of cytokine-mediated seeding and marginal chemotaxis (Fig. 8). Seeding cells at the currently available maximum (i.e., a sampling fraction of 1) produced the least effective immune response to prevent infection for a given chemotaxis parameter, as was marginal chemotaxis for a given sampling fraction. Likewise, seeding immune cells mostly at random produced a better strategy for mitigating spread of infection as did stronger chemotactic sensitivity for a given sampling fraction. Immune response was most effective for sampling fractions up to 0.1% and was more effective with increasing sampling fraction and decreasing chemotaxis parameter.

Fig. 8
figure 8

Susceptible cells from 100 replicas of the spatial model of viral infection while varying immune cell sampling fraction and chemotaxis parameter. Medians shown as black lines. 0th to 100th quantiles shaded as blue, 10th to 90th as orange, and 25th to 75th as light blue. ODE model results shown as red dashed lines. Star marks best. Triangle marks worst. Circles with green to red shading show best to worst total error, respectively

A significant source of differences between the spatial and ODE models was differences in the timing of peak viral load during progression of infection (Fig. 9). For greater sampling fractions and lesser chemotaxis parameters, peak viral load was greater and occurred later in replicas using the spatial model. These differences produced viral load curves for simulation replicas of the spatial model with the same dynamical features as the ODE model (maximum value around 2 to 3 days, followed by logarithmic decay to zero), and even (nearly) identical results before peak viral load, but with continued progression of infection up to a day longer (e.g., sampling fraction of 1.0 and chemotaxis parameter of 104, worst and best QRF 12% and 60%, dead cells and infected cells, respectively; worst and best NMAE 53% and 17%, susceptible cells and local immune cells, respectively) before decline and elimination of infection. For mostly random seeding and strong chemotaxis (e.g., sampling fraction of 2.5 × 10−5 and chemotaxis parameter of 106), the only notable differences between spatial and ODE model results occurred long after peak viral load and for all parameter sets, where single infection events produced strong fluctuations in otherwise negligible viral load values. The overall best-fit parameter set was a sampling fraction of 0.01% and chemotaxis parameter of 106 (worst and best QRF 24.5% and 82%, lymph immune cells and infected cells, respectively; worst and best NMAE 24% and 5%, virus-releasing cells and dead cells, respectively).

Fig. 9
figure 9

Viral load from 100 replicas of the spatial model of viral infection while varying immune cell sampling fraction and chemotaxis parameter. Medians shown as black lines. 0th to 100th quantiles shaded as blue, 10th to 90th as orange, and 25th to 75th as light blue. ODE model results shown as red dashed lines. Star marks best. Triangle marks worst. Circles with green to red shading show best to worst total error, respectively

Lastly, we tested the best-fit parameters for aforementioned parameter sweep for the scenario of one initially infected cell. Previously observed stochasticity of system dynamics and outcomes for one initially infected cells in “Two-dimensional infectivity” were also observed when using the best-fit parameters (Fig. 10). In the case of also simulating an immune response, the final distribution of local immune cells varied significantly among 100 simulation replicas. In simulation replicas that experienced more infection, few immune cells were found near the site of initial infection, whereas in replicas that experienced less infection, more immune cells were found nearer to the site of initial infection.

Fig. 10
figure 10

Spatial model results from two simulation replicas of viral infection and immune response in a quasi-two-dimensional, epithelial sheet with best-fit immune response parameters and one initially infected cell. Results shown at 0, 3, 4, 5, 7 and 14 days in simulation time. Epithelial and immune cells shown as in Fig. 6

Two simulation replicas experienced no significant infection, and those replicas that did experienced infection at significantly varying degrees of severity (Fig. 11). Among replicas with more than one infected cell during simulation time, the final number of susceptible cells was in the range of 772 to 1045 cells, with mean and standard deviation of 880 and 98.7 cells, respectively. Compared to a final number of susceptible cells of 722.7 according to the ODE model, the spatial model produced results that disagree with the ODE model in the range of about 7 to 45% less infection for replicas that experienced infection (worst and best QRF 29% and 82%, susceptible cells and virus-releasing cells, respectively; worst and best NMAE 91% and 16%, local immune cells and susceptible cells, respectively). Measurements of viral load and local and lymph cytokine were comparable to the ODE model in magnitude, though the replicas that experienced less severe infection also experienced later peak viral loads and cytokine, with corresponding effects on immune cell recruitment.

Fig. 11
figure 11

Scalar results from 100 replicas of the spatial model of viral infection and immune response with best-fit immune response parameters and one initially infected cell. Medians shown as black lines. 0th to 100th quantiles shaded as blue, 10th to 90th as orange, and 25th to 75th as light blue. ODE model results shown as red dashed lines

Discussion

We demonstrated that spatial mechanisms like mass transport are implicitly represented in non-spatial model parameters like virus infection rate. Referring to Eq. 1, the model term (−βUV) describes the rate at which susceptible cells become infected by exposure to virus without regard for which susceptible cell is infected, or for which virion does the infecting, or for how the virion arrived at the cell. The ODE model uniformly exposes each cell to each virion and employs the model infectivity parameter β to describe all other mechanisms associated with infection besides the presence of susceptible cells and virions, including transport of virions from virus-releasing cells to susceptible cells. When varying the diffusion coefficient of virus, we found that we could map reductions in virus diffusion to reductions in the ODE model infectivity, which produced ODE model results that were indistinguishable from spatial model median results for all but the most extreme parameter sets (i.e., very low diffusion and few initially infected cells, Figs. 4 and 5). Referring to our named outcomes in “Conclusion.”

In this work, we developed a method for generating spatial, multicellular models of biological systems from non-spatial models, and vice versa, which we call cellularization. We demonstrate using our method by cellularizing non-spatial models of viral infection and host-pathogen interaction. Using these cellularized models, we quantitatively showed how spatial mechanisms implicitly represented in non-spatial models can exhibit significant effects on emergent dynamics when explicitly modeled. Variations in related non-spatial model parameters emerged from moderate cases of varying spatial mechanisms like rate of diffusive mass transport of virus, while extreme cases generated emergent dynamics inconsistently with those described by the non-spatial model. We describe the responsible mechanisms for extreme disagreement between homogeneous and cellular models, specifically concerning the limitations of describing discrete biological objects and processes using continuous descriptions.

Methods, we found that sufficiently high virus diffusion results in an Ensemble average, sufficiently low virus diffusion results in incompatibility, and between these two extrema results in localization. This aspect of a non-spatial model is particularly important when considering how to interpret in vitro results and translate information gained from in vitro scenarios to predictions of in vivo outcomes. For example, diffusive transport of extracellular virions during viral infection depends on the medium through which virions migrate, which varies in vitro by culture medium and in vivo by location. Results in “Effects of diffusivity” can be regarded in this context as characterizing the effects of the environment on progression of infection while holding all other mechanisms constant, of which we showed to significantly affect both the rate and severity. Such observations then elucidate a means by which we can interrogate when coarse-grained approaches fail to reliably predict physical data by examining the underlying spatiotemporal mechanisms. As demonstrated in this work, our method of cellularization enables such activities related to both synthesis and validation of non-spatial models. Furthermore, the model parameters for the diffusion coefficient of cytokines and virus in immune response simulations were of comparable magnitude, though cytokines are often reported to be an order of magnitude greater than that of virus. While the purpose of this work is not to generate accurate predictions of a particular virus and cytokine in a particular location, it is interesting to note that the reduction in magnitude of diffusion of small molecules alone has been observed in human cervical mucus, where larger viruses diffuse unhindered [30, 31].

We also demonstrated how a generated spatial analogue of a non-spatial ODE model using cellularization can further elucidate the features of the underlying mechanisms that a non-spatial model only implicitly represents. The parameter sweep of spatial model parameters in “Immune response modeling” showed the possible strategy of the immune response to eliminate viral infection. Specifically, immune cells more effectively induce apoptosis in infected cells via contact-mediated interactions when they arrive at a mostly random location near a local site of infection and then strongly chemotax along soluble immune signals. Probably, arriving where immune signals are the strongest results in immune cells that cannot effectively prevent the advancing front of an infection, due to the nature of diffusive transport. Rather, it may be more effective to arrive near, but just outside, a lesion. This suggests that resident cells of a site of infection may further activate chemotactic sensitivity in the immune cells described by an ODE model like ours that models viral infection and immune response (e.g., CD8+ T cells). Such observations well describe recruitment and search strategies within the context of current understanding of effector T cells [32]. Namely, that effector T cell search strategies for infected cells balance a trade-off of exploration and exploitation (i.e., of random and directional migration), and that T cells presumably encounter antigen-presenting cells much more often in inflamed tissue than in lymph nodes, which changes their search strategy from exploration-dominant to exploitation-dominant. Cellularization of detailed ODE models of the immune response [25] could generate detailed multiscale, spatial, cell-based models, which would otherwise be especially difficult to construct de novo, that include these details for the quantitative study of the immune response with spatial information and cellular precision.

Some simulation replicas of viral infection generated no infection of tissue even without an immune response for various model parameters and initial conditions. This observation is particularly important with regard to both the biological features represented in the non-spatial and spatial models, as well as the meaning of spatial model, of viral infection. It is certainly a reasonable and relevant question to ask whether a single infected cell can cause infection in tissue, and to what extent. For such a question, the ODE model always predicts infection, due to its continuous treatment of discrete objects and processes (e.g., a half of a virion can generate a quarter of an infected cell in the ODE model), while the spatial analogue allows for premature death and ineffective viral exposure. In this regard, the spatial analogue better describes the total story of infection. This is not to say that either is a sufficiently insightful or predictive model of viral infection per se, or that failure to infect is unique to the generated spatial model in general [17, 33]. Rather, we argue that the spatial analogue provides a better description of the biophysical mechanisms responsible for viral infection according to the ODE model. The feasibility of employing such a spatial analogue on useful modeling applications could be questioned due to the greater computational cost of multicellular simulations, though such a limitation is technological. To this end, our cellularization method improves the feasibility of using spatial, multicellular models by readily permitting the employment of model parameters fitted to in vivo and/or in vitro data using analogous non-spatial models, which are computationally inexpensive.

More importantly, deciding on the context of the spatial domain of a spatial analogue with respect to the entire biological system is a particularly challenging modeling problem. For example, the spatial analogues developed in this work employed periodic boundary conditions on the boundaries orthogonal to the epithelial sheet. This arrangement implies that the spatial analogue describes one instance of a periodic system, the collection of which the ODE model describes. This paradigm breaks down when considering highly localized cases like one initially infected cell, as demonstrated by multiple simulation replicas of the spatial analogue producing significantly different outcomes (i.e., widespread infection or no infection). In such a case, we could observe no infection in one replica and then argue that infection occurred in some other replica, at the expense of revising the premise of the model altogether (i.e., the system described by the ODE model is not periodic). As inconvenient as such a scenario may be, it demonstrates that cellularization better equips the modeler to interrogate the underlying mechanisms responsible for the emergent dynamics of a biological system by providing a means by which explicit cell-based and/or spatial models can be consistently generated from non-spatial models. In the case of viral infection, cellularization enables the utilization of both well-fitted ODE models of in vivo data and in vitro evidence of viral transport and discrete viral internalization by inquiring about the cellular, multicellular, and spatiotemporal mechanisms that can consistently explain observations at multiple scales. This approach is particularly powerful because it enables mechanistic modeling directly within the context of emergent dynamics and on the basis of individual cells, the latter of which is neglected in common, homogenized spatial modeling approaches like reaction-diffusion and traveling-wave models. Whether the mathematical formalism of cellularization can be transformed into a compatible structure with such modeling approaches is currently unclear, though the ability to translate models between those approaches and cell-based approaches clearly be valuable to biological modeling.

Future work

In the immune response model, the spatial model produced greater lymph immune cells and fewer local immune cells than the ODE model. These differences may be due in part to that the seeding algorithm always considers increments in the number of cells under consideration for seeding beginning at zero. The significance of this apparent bias is currently unclear, as is any exact mitigation strategy, though it may be possible to refine the seeding algorithm by adding stochasticity to the considered number of seeded cells per algorithmic iteration.

There are many mathematical terms employed by non-spatial models of biological systems to which the general forms described in this work (e.g., Eqs. 21, 26 and 30) cannot be readily applied to generate linear stochastic transition rules during cellularization (e.g., Michaelis-Menten). Deriving spatial analogues of ODE model terms not only provides broader utility by allowing the development of an analogous spatial model, but also a stronger understanding of the mechanism(s) represented by the ODE model term. Likewise, proving the nonexistence of a spatial analogue calls into question what exactly an ODE model term means. In such a case for an ODE model term that describes a cellular interaction, questions arise as to what a cellular model means if it cannot be written on a cellular basis. This is not to propose that the validity of an ODE model depends on its mathematical compatibility with cellularization. Rather, we propose that an ODE model term that is incompatible with cellularization represents multiple mechanisms, each of which can be described by an ODE model term that supports cellularization. The justification for this proposition is obvious from a biological perspective, in that any subsystem of an organism defined at a scale higher than the cellular scale can be resolved to a cellular basis, which is the basis of spatial models according to our cellularization method. A clear separation of mechanisms provides further insights into the dynamics of a biological system relevant to both biomedical sciences (e.g., when developing drug therapies) and basic biology (e.g., when quantitating biological complexity), as well as a more descriptive and robust mathematical model.

We can also turn this paradigm—of regarding ODE models within the context of generating a spatial analogue—towards spatial models of multicellular systems while considering their non-spatial analogue. Furthermore, operations defined in our cellularization method on spatial models may generate novel non-spatial model terms that are not obvious from a homogenized approach. This approach may also be useful for introducing better descriptions of such spatial effects as those observed in “Effects of diffusivity” to non-spatial models, or for estimating non-spatial model parameters from in vivo imaging data in clinical applications. To this end, cellularization provides the ability to generate non-spatial models to describe various “compartments” of a biological system (e.g., upper and lower respiratory tract compartments), where the non-spatial model well describes the underlying biocomplexity of a compartment according to its local specificity. We envision broad application of cellularization to do detailed spatiotemporal modeling of a particular biological module and then propagate information to, from, and across coarser scales. Indeed, previous work has cast such a vision concerning the hierarchical organization and function of the liver [34]. Cellularization provides a clear and consistent approach to develop such multiscale modeling frameworks of whole organisms with regard to the scale of the cell.

Conclusion

In this work, we developed a method for generating spatial, multicellular models of biological systems from non-spatial models, and vice versa, which we call cellularization. We demonstrate using our method by cellularizing non-spatial models of viral infection and host-pathogen interaction. Using these cellularized models, we quantitatively showed how spatial mechanisms implicitly represented in non-spatial models can exhibit significant effects on emergent dynamics when explicitly modeled. Variations in related non-spatial model parameters emerged from moderate cases of varying spatial mechanisms like rate of diffusive mass transport of virus, while extreme cases generated emergent dynamics inconsistently with those described by the non-spatial model. We describe the responsible mechanisms for extreme disagreement between homogeneous and cellular models, specifically concerning the limitations of describing discrete biological objects and processes using continuous descriptions.

Methods

In this section, we describe our method for converting model objects and processes between non-spatial ODE models and spatial, cell-based models of multicellular systems. We consider paired models of the same underlying biological/biochemical system, where an ODE model consists of only scalar-valued variables related by rate laws that depend only on the variables and a set of scalar parameters, and hence is intrinsically non-spatial. The spatial model could include scalar variables, continuously variable spatial quantities (i.e., fields) and discrete agents with individual states, which typically have spatial locations (which could change in type) and which can come into existence (e.g., birth) and disappear (e.g., death). The relationships among agents and variables are more diverse in spatial models, including the rate equations of the non-spatial case, as well as stochastic transition rules, partial-differential equations describing the effect of spatial variation on rates, and more complex boundary and initial condition considerations. Initial conditions require consistent specification of initial values of variables only among paired models, with a limited number of possible boundary conditions in the spatial model. Rate laws with explicit time delays and integro-differential forms are extensions of the basic forms treated in this work, though many of the issues we consider apply. Likewise, non-spatial agent-based models form an intermediate class, which we do not consider in detail here, though many of the same considerations apply.

We define a spatial model as analogous to a non-spatial ODE model if the ODE model describes the dynamics of the spatial model in the limit of well-mixed conditions. As such, we refer to a spatial model generated from an ODE model as the spatial analogue of the ODE model. In this work, we consider cells, diffusive species, and the interactions between them. We assume that both an ODE model and its analogous spatial model describe the same underlying real-world dynamics of a biological system, but that only the spatial model explicitly describes spatial dynamics. We assume that the ODE model assumed well-mixed conditions such that all ODE model variables, parameters, and equations represent spatially homogeneous properties of, and interactions between, underlying objects. Likewise, we assume that the spatial domain in which the objects and processes are described by an analogous spatial model represents a constituent spatial element of the biological system described by the ODE model. Under these assumptions, the following outcomes are possible when comparing results from a pair of ODE and spatial models:

  • Homogeneity. The models agree for the same parameters and the system is spatially homogeneous.

  • Ensemble average. The models agree for the same parameters and the system is spatially heterogeneous.

  • Localization. The models agree for different parameters and the system is spatially heterogeneous.

  • Incompatibility. The models never agree.

We define two metrics to quantify the degree of similarity between the ODE model time series and the ensemble of cellularized model time series for any scalar quantity predicted by both models. The first metric, the normalized mean absolute error (NMAE), is the time average of the absolute value of the difference between the ODE model prediction at each time and the median over all cellularized model replica predictions at that time point. This metric would be zero if the ODE and median cellularized model time series agreed exactly (identical values at each time) and increases in an intuitive way with visual distance between the two curves. The second metric, the quantile range fraction (QRF), is the fraction of the time during which the value predicted by the ODE model at a given time falls within a given range of quantiles for the ensemble of cellularized model prediction. We use the range of 10th to 90th quantiles in our QRF. The QRF quantifies the degree to which the ODE model time series falls within the distribution of time series for the ensemble of cellularized model replicas. The second metric is less sensitive to specific values of the quantity examined than the NMAE.

While in general nothing prohibits us from defining multiple spatial domains that represent various interesting elements of the biological system, we assume the simplest case in this work when deriving analogous mathematical forms—that a spatial domain is a representative volume element of the entire biological system described by an ODE model. We describe our formalism as beginning with an existing ODE model and deriving an analogous spatial model. However, our formalism works in the opposite direction, allowing us to derive an ODE model from an existing spatial model.

Spatiotemporal scaling using well-mixed conditions

Well-mixed conditions allow us to relate measurements of quantities at different scales, where measurements at every scale everywhere exactly and uniformly scale. For a measurement Z of a quantity (e.g., number of cells of a particular type, total amount of a species in a volume) at one scale and z measuring the quantity of the same object but at another scale, let μ be a coefficient relating the two scales such that

$$ z=\mu Z. $$
(7)

The same is true for the behavior of objects according to well-mixed conditions, where changing the scale at which an object operates does not qualitatively affect the dynamics of the object described by the ODE model. As such, for a rate equation of the form,

$$ \frac{dZ}{dt}=F\left(Y,Z\right), $$
(8)

dynamics for a measurement z at some other scale and y = μY take the form,

$$ \frac{dz}{dt}\left(Y,Z;\mu \right)=\mu F\left(\frac{y}{\mu },\frac{z}{\mu}\right). $$
(9)

Scaling of ODE model quantities under well-mixed conditions is assumed to be true for quantities of mass and number (e.g., the population of cells of a particular type) at every scale. Depending on the scale of a multiscale spatial model, well-mixed conditions allow us to cast a quantity Z from an ODE model as acting globally (i.e., as spatially homogeneous) or locally (i.e., as spatially heterogeneous) in a spatial model, depending on the spatial qualities imposed upon Z. For example, diffusive species with very high or very low diffusivity may be cast as acting globally or locally, respectively, at the scale of individual cells. Henceforth, for a measurement Z of the quantity of an object at the scale of an ODE model, let z be a global (i.e., non-spatial) measurement of the quantity of the same object at the scale of a spatial model, and let \( \overset{\sim }{z} \) be an infinitesimal, point-like measurement of the same object.

To employ Eq. 9 to describe the dynamics of a quantity z at the scale of a spatial domain according to an ODE model requires the relating of a static measurement \( {\mathcal{M}}_o \) of some scalar quantity at the scale of the ODE model to a global, static measurement \( {\mathcal{M}}_s \) of a scalar quantity at the scale of its analogous spatial model such that, using Eq. 7, for a global scaling coefficient η and time t,

$$ \eta :{\mathcal{M}}_o\mapsto {\mathcal{M}}_s\forall t. $$
(10)

For example, an ODE model describing a volume \( \left\Vert {\mathcal{U}}_o\right\Vert \) and an analogous spatial model describing a volume \( \left\Vert {\mathcal{U}}_s\right\Vert \) presents a form of Eq. 10 with a scaling coefficient ηvol = η,

$$ \left\Vert {\mathcal{U}}_s\right\Vert ={\eta}_{vol}\left\Vert {\mathcal{U}}_o\right\Vert . $$
(11)

Likewise, an ODE model describing some fixed number \( {\mathcal{T}}_o \) of a set of cell types and analogous spatial model describing a fixed number \( {\mathcal{T}}_s \) of a set of the same cell types presents a form of Eq. 10 with a scaling coefficient ηcount = η,

$$ {\mathcal{T}}_s={\eta}_{count}{\mathcal{T}}_o. $$
(12)

In the case of Eq. 11, application of Eq. 9 to generate global measurements at the scale of a spatial model would presuppose that, under well-mixed conditions, a quantity Z is the same per volume. Likewise for Eq. 12, well-mixed conditions would impose that a quantity Z is the same per cell.

For the case of counting numbers of a cell type in a discrete cell population, let σ be the identifier of a cell and let τ(σ) be the type of σ. Let an ODE model of the quantity of a cell type \( \hat{N} \) be written as N. Using this convention, the discrete analogue of N naturally follows by counting the number of cells of type \( \hat{N} \),

$$ N=\left\Vert \left\{s:\tau (s)=\hat{N}\right\}\right\Vert . $$
(13)

Likewise, for local measurements of some global, scalar measurement Z at some scale, consider a local, point-like measurement \( \overset{\sim }{z}=\overset{\sim }{z}\left(x,t\right) \) of Z in a space \( \mathcal{U}\subset {\mathbb{R}}^n \), where \( x\in \mathcal{U} \). For diffusive \( \overset{\sim }{z} \), let \( \overset{\sim }{z} \) be governed by a reaction-diffusion equation of the general form,

$$ {\partial}_t\overset{\sim }{z}={\partial}_i\left({D}_Z{\partial}_i\overset{\sim }{z}\right)+\overset{\sim }{v}\overset{\sim }{z}+\overset{\sim }{w}. $$
(14)

Considering the case of homogeneous \( \overset{\sim }{z} \) (i.e., under well-mixed conditions), Z and \( \overset{\sim }{z} \) can be related with a local scaling coefficient θ using Eq. 7,

$$ \underset{D_Z\to \infty }{\lim}\overset{\sim }{z}=\theta Z, $$
(15)

and in general, for global quantity z at the scale of a spatial domain (i.e., z = ηZ),

$$ z={\int}_{\mathcal{U}}\overset{\sim }{z}(x) dV. $$
(16)

Then, when considering Eq. 16 for a homogeneous \( \overset{\sim }{z} \),

$$ \theta =\frac{\eta }{\left\Vert \mathcal{U}\right\Vert }. $$
(17)

When deriving a spatial model from ODEs, the following relationships are assumed to hold under well-mixed conditions for ODE model measurement Z, spatial global measurement z, and local, point-like measurement \( \overset{\sim }{z} \), according to Eqs. 11 and 17

$$ Z=\frac{1}{\eta }z=\frac{1}{\theta}\underset{D_Z\to \infty }{\lim}\overset{\sim }{z}.. $$
(18)

Using Eq. 18, quantities described by an ODE model at some scale can then be employed in a spatial model at some smaller scale as either globally acting objects (according to the ODE model, using Eq. 9) or a locally acting object (according to Eq. 14 for diffusive species, and to subsequent derivations for cells).

To define a point-like measurement of a quantity associated with a particular cell, well-mixed conditions are also applied to the domain of a cell. Let σ(x, t) be a cell identifier that describes the location of individual cells in the spatial domain (e.g., for cell s, σ(x, t) = s at every point x occupied by cell s at time t), and let \( \mathcal{V}\left(s,t\right)=\left\{x\in \mathcal{U}:\sigma \left(x,t\right)=s\right\} \) be the domain of cell s at time t in the spatial domain \( \mathcal{U}\subset {\mathbb{R}}^n \). A point-like measurement \( \overline{z}\left(s,t\right) \) of a spatially heterogeneous quantity \( \overset{\sim }{z} \) associated with a cell s at time t under well-mixed conditions is then

$$ \left\Vert \mathcal{V}\left(s,t\right)\right\Vert \overline{z}\left(s,t\right)={\int}_{\mathcal{V}\left(s,t\right)}\overset{\sim }{z}\left(x,t\right) dV. $$
(19)

Note that for point-like cellular objects,

$$ \underset{\left\Vert \mathcal{V}\left(s,t\right)\right\Vert \to 0}{\lim}\overline{z}\left(s,t\right)=\overset{\sim }{z}\left(x,t\right)\forall x\in \mathcal{V}\left(s,t\right). $$
(20)

Cellularization of signals

Consider a rate equation for chemical species F and G and number N of a cell type \( \hat{N} \),

$$ \frac{dF}{dt}=v\left(F,G\right)+w\left(F,G\right)N. $$
(21)

Here v and w are generic functions with arguments F and G for demonstrative purposes, but are otherwise arbitrarily selected (i.e., they could be any set of chemical species or some other object of the system). Using Eq. 18, let \( \overset{\sim }{f} \) be the heterogeneous analogue of F according to Eq. 14 with diffusivity DF. In general, v can be applied in any spatially particular way (e.g., only as a flux on a boundary) on the condition that the volume integral of its spatial analogue is equal to ηv at every time t. We consider here the case when v is uniformly applied in the spatial domain,

$$ {\partial}_t\overset{\sim }{f}\left(x,t\right)={\partial}_i\left({D}_F{\partial}_i\overset{\sim }{f}\left(x,t\right)\right)+\theta v\left(\frac{\overset{\sim }{f}\left(x,t\right)}{\theta },\frac{\overset{\sim }{g}\left(x,t\right)}{\theta}\right)+\overset{\sim }{w}\left(x,t\right). $$
(22)

Here v is rewritten as a substitution of \( \overset{\sim }{v}\overset{\sim }{z} \) in Eq. 14 using Eq. 9. In the case of heterogeneous \( \overset{\sim }{f} \), \( \overset{\sim }{w} \) can be related to wN in Eq. 21 by integrating over all cell volumes of type \( \hat{N} \),

$$ w\left(F(t),G(t)\right)N(t)=\sum \limits_{s\in \left\{{s}^{\prime }:\tau \left({s}^{\prime },t\right)=\hat{N}\right\}}{\int}_{\mathcal{V}\left(s,t\right)}\overset{\sim }{w}\left(x,t\right) dV. $$
(23)

We can write an expression for \( \overset{\sim }{w} \) at each point as a function of \( \overset{\sim }{f} \) and \( \overset{\sim }{g} \) such that Eq. 23 is true under well-mixed conditions by employing Eq. 19,

$$ \overset{\sim }{w}\left(x,t\right)=\left\{\begin{array}{cc}\frac{1}{\left\Vert \mathcal{V}\left(\sigma \left(x,t\right),t\right)\right\Vert }w\left(\frac{\overline{f}\left(\sigma \left(x,t\right),t\right)}{\theta },\frac{\overline{g}\left(\sigma \left(x,t\right),t\right)}{\theta}\right)& \tau \left(\sigma \left(x,t\right),t\right)=\hat{N}\\ {}0& \mathrm{otherwise}\end{array}\right.. $$
(24)

Note that w in Eq. 19 and \( \overset{\sim }{w} \) in Eq. 24 are general in the sense that multiple terms could be defined with respect to multiple cell types. In such a case, each type-specific term would have a corresponding \( \overset{\sim }{w} \) according to the general form in Eq. 23.

Hybridization of population dynamics

We derive discrete, stochastic events from ODE rate equations by considering the ODE rate equations as mean-field approximations of discrete, stochastic events. The total number of occurrences of a discrete, stochastic event K is assumed to be a random variable drawn from a Poisson distribution, which allows us to relate expressions from ODE and discrete, stochastic models. We use the cumulative distribution function \( {\mathcal{P}}_k\left(r\Delta t\right) \) for determining the probability that a discrete event with mean rate r occurs more than k times over a simulation step of period ∆t,

$$ {\mathcal{P}}_k\left(r\Delta t\right)=\Pr \left(K>k\right)=1-{e}^{-r\Delta t}\sum \limits_{0\le n\le k}\frac{{\left(r\Delta t\right)}^n}{n!}. $$
(25)

For a discrete event that can occur exactly once in a simulation step (e.g., killing of one cell), we use \( {\mathcal{P}}_0 \), which determines the probability of an event having occurred at least once (rather than exactly once).

Consider ODE model population dynamics of the number of cell types \( \hat{M} \) and \( \hat{N} \) of the form,

$$ \left\{\begin{array}{c}\frac{dN}{dt}=f- gN- uN\\ {}\frac{dM}{dt}= uN\end{array}\right.. $$
(26)

Three events are described in the rate equation for N. First, the term f describes the net inflow of \( \hat{N} \)-type cells. Second, the term −gN describes the net outflow of \( \hat{N} \)-type cells. Third, the term −uN describes the transition of a cell from a \( \hat{N} \)-type cell to a \( \hat{M} \)-type cell (e.g., from a living cell to a dead cell), which we denote as \( \hat{N}\to \hat{M} \). Discrete, stochastic analogues are derived in this order using a mean-field approximation of the Poisson distribution.

The net inflow of \( \hat{N} \)-type cells is written using the cumulative distribution function of the Poisson distribution for the occurrence of k cells of type \( \hat{N} \) being added to the spatial domain,

$$ \Pr \left(\mathrm{add}\ k\ \hat{N}-\mathrm{type}\ \mathrm{cells}\right)={\mathcal{P}}_k\left(f\Delta t\right). $$
(27)

Eq. 27 is implemented using the following algorithm. Beginning with k = 0, draw a uniformly distributed random number X in [0, 1]. If X is greater than Eq. 27 for k, then add k cells of type \( \hat{N} \). Otherwise, increment k by one and repeat. Functionally, this imposes that, while counting upwards from k = 0, if no more than k cells are added, then k cells are added.

The net outflow of \( \hat{N} \)-type cells is considered for each \( \hat{N} \)-type cell independently at each simulation step as an event that can occur no more than once. For cell s,

$$ \Pr \left(\mathrm{remove}\ s|\tau (s)=\hat{N}\right)={\mathcal{P}}_0\left(g\Delta t\right). $$
(28)

Note that a form like Eq. 28 can also be generated for describing the probability of mitosis in the case where the term −gN in Eq. 26 were instead +gN.

The discrete transition event of a cell of type \( \hat{N} \) to type \( \hat{M} \) over a period ∆t is considered for each \( \hat{N} \)-type cell of identification s and type τ(s, t) at time t as an event that can occur no more than once,

$$ \Pr \left(\tau \left(s,t+\Delta t\right)=\hat{M}|\tau \left(s,t\right)=\hat{N}\right)={\mathcal{P}}_0\left(u\Delta t\right). $$
(29)

Hybridization of contact-mediated population interactions

Contact-mediated interactions are critical to many developmental (e.g., stem cell proliferation) and physiological (e.g., antigen presentation) processes. In addition to the interactions described in “Hybridization of population dynamics,” suppose that u in Eq. 26 describes a rate of \( \hat{N}\to \hat{M} \) mediated by contact interactions between \( \hat{N} \)- and \( \hat{O} \)-type cells and, for a > 0 and cell types \( \hat{M} \), \( \hat{N} \), and \( \hat{O} \), Eq. 26 has the form,

$$ \left\{\begin{array}{c}\frac{dM}{dt}= aNO+f\\ {}\frac{dN}{dt}=- aNO+g\end{array}\right.. $$
(30)

Spatial models of contact-mediated population interactions are derived by considering well-mixed conditions in the strictest sense. If \( \hat{M} \)- and \( \hat{N} \)-type cells constitute a cellular body and present a certain surface area available for contact interfaces, then the surface area of contact interfaces between \( \hat{N} \)- and \( \hat{O} \)-type cells in well-mixed conditions is proportional to the number of \( \hat{N} \)- and \( \hat{O} \)-type cells N and O, respectively.

As such, let AP be the total available contact surface area of all \( \hat{P} \)-type cells of fixed number, let P = M + N be a quantity of a set of \( \hat{M} \) and \( \hat{N} \) cells, let AN be the total available contact surface area of a \( \hat{N} \)-type cell, and let AN, O be the contact area of a \( \hat{N}-\hat{O} \) interface. Under well-mixed conditions,

$$ {A}_{N,O}=\frac{O{A}_O}{A_P}{A}_N. $$
(31)

Eq. 31 can be rewritten for O, in which case Eq. 30 can be rewritten on a cellular basis for cell s with total interfacial contact area with \( \hat{O} \)-type cells As, O (i.e., AN, O → As, O on a cellular basis). \( \hat{M}\to \hat{N} \) mediated by \( \hat{N}-\hat{O} \) contact interfaces is then considered for each cell of identification s and type \( \hat{N} \) as an event that can occur no more than once,

$$ \Pr \left(\tau \left(s,t+\Delta t\right)=\hat{M}|\tau \left(s,t\right)=\hat{N}\right)={\mathcal{P}}_0\left(\gamma (s)\Delta t\right), $$
(32)

where

$$ \gamma (s)=\frac{a{A}_P}{A_N}\frac{A_{s,O}}{A_O}. $$
(33)

Eq. 33 scales a from the general ODE form in Eq. 30 to the ratio of two ratios, namely the ratio of the contact area of a cell with \( \hat{O} \)-type cells to the total contact area of a \( \hat{O} \)-type cell, to the ratio of the total contact area of a \( \hat{N} \)-type cell to the total contact area of all \( \hat{P} \)-type cells. Note that in the case of geometrically identical \( \hat{P} \)-type cells (i.e., \( \hat{M} \)- and \( \hat{N} \)-type cells have the same total contact surface area), Eq. 33 can be rewritten to scale a to the total number of \( \hat{P} \)-type cells,

$$ \gamma (s)= aP\frac{A_{s,o}}{A_o}\leftrightarrow {A}_P={PA}_N $$
(34)

Hybridization of recruitment

Eq. 26 describes the inflow of cells with a rate f and outflow of cells with a rate g per cell. Here we describe handling inflow of cells, which can be generally described, and neglect describing any general method for outflow, which is particular to problem setup and cellular dynamics method (e.g., whether cells are simply removed wherever they are when outflow occurs, or they first migrate towards a boundary). To describe the inflow of cells in a spatial model, it is necessary to define a boundary through which each incoming cell enters the spatial domain, and the location on that boundary where each incoming cell appears. For an inflow rate f into a spatial domain \( \mathcal{U} \), we refer to a corresponding boundary \( {\mathcal{W}}_f\subseteq \partial \mathcal{U} \) on which cells appear as an inflow boundary. Likewise, we refer to placing a new cell on an inflow boundary as seeding. A cell can be seeded at a site on an inflow boundary that is not occupied by another cell. We refer to such sites as being available and denote the set of available sites on an inflow boundary for an inflow rate f as \( \partial {\mathcal{W}}_{f,0} \),

$$ \partial {\mathcal{W}}_{f,0}(t)=\left\{x\in \partial {\mathcal{W}}_f:\sigma \left(x,t\right)=0\right\}. $$
(35)

In general, we define an available site selection function \( \mathcal{S}\left(\zeta, t\right) \) that selects an available site for seeding a cell of type ζ at time t. It should be noted that an available site selection function does not describe the probability of seeding a cell (e.g., Eq. 27), but rather maps the set of available sites on inflow boundaries to a coordinate where seeding a cell occurs for a given cell type and time. It should also be noted that an available site selection function is defined to take the argument of a cell type, rather than an inflow rate, since one could define a total inflow rate as consisting of multiple inflow rates on multiple inflow boundaries (e.g., two thirds of incoming cells enter on one boundary, and one third of incoming cells enter on another boundary). Imposing no additional spatial information on seeding an incoming cell results from designating all boundaries of the spatial domain as inflow boundaries, and then randomly selecting an available site with a uniform distribution. Likewise designating a subset of the boundaries of a spatial domain as inflow boundaries necessarily imposes additional spatial information onto the spatial analogue of an ODE model by introducing spatial heterogeneity to the boundaries of the spatial domain.

An available site selection function of particular interest (though not critical to the total method presented here) that we propose concerns the influence of local heterogeneity of signals in a spatial domain that induce directional migration (e.g., chemotaxis, haptotaxis) in incoming cells. Suppose that a locally heterogeneous field c(x, t) attracts an incoming cell of type ζ that is to be seeded somewhere in a set of available sites \( {\partial \mathcal{W}}_{g,0}(t) \) on an inflow boundary according to a rate g. Rather than impose additional spatial information outside of a spatial domain, we can instead presuppose that c near \( {\partial \mathcal{W}}_{g,0}(t) \) influences an otherwise random seeding location where an incoming cell of type ζ arrives. In this case, consider a randomly selected subset of available sites \( {\omega}_{g,0}(t)\subseteq {\partial \mathcal{W}}_{g,0}(t) \), where \( \left\Vert {\omega}_{g,0}(t)\right\Vert ={w}_g\left\Vert {\partial \mathcal{W}}_{g,0}(t)\right\Vert \) for a sampling fraction wg (0, 1] and each selected site yωg, 0(t) is selected with a uniform probability (i.e., without imposing any additional spatial information, as mentioned in the preceding text). We impose the effects of a locally heterogeneous attractant c of cell type ζ at time t by seeding an incoming cell of type ζ at the location where c is greatest among the randomly selected subset of available sites,

$$ \mathcal{S}\left(\zeta, t\right)\in \underset{y\in {\omega}_{g,0}(t)}{\mathrm{argmax}}c\left(y,t\right). $$
(36)

It naturally follows that replacing argmax(∙) with argmin(∙) imposes a repellant c. Note that additional selections must be made when this particular available site selection function yields multiple values (e.g., we randomly choose a site when c is everywhere zero). Also note that wg = 1 imposes that an incoming cell is seeded exactly at the location with the maximum value of c in all available sites \( {\partial \mathcal{W}}_{g,0}(t) \), while wg → 0 imposes that an incoming cell is seeded at a random location.

Mixing mixed conditions

Nothing prohibits us from supposing that only some of the cells of a population act locally, while the remaining cells of the population act globally, in the sense that only the locally acting cells are represented in a spatial domain, while the effects of the entire population are present in the spatial domain. We can, for example, say that of 2N cells of type \( \hat{N} \), only N cells are present in a spatial domain, the interactions of which are modeled locally, while the other N cells are not present in the spatial domain, and their interactions occur everywhere homogeneously (in fact, the well-mixed conditions employed by ODE models necessitate that this is true). The same is true concerning signal fields, in that a homogeneous field can either be explicitly represented in the spatial domain, or represented as a scalar value equal to the volume integral of an equivalent homogeneous field. This distinction is important for at least two reasons.

First, from a computational standpoint, explicitly integrating the diffusion equation in time becomes increasingly computationally expensive as the diffusion coefficient (e.g., DF in Eq. 22) of a soluble signal increases (because of numerical stability). Likewise, as the diffusion coefficient of a soluble signal increases, the field acts more homogeneously. Other factors permitting (e.g., diffusion length), one can mitigate the computational cost of integrating fast-acting soluble signals in time by representing them as scalar-valued functions (according to their ODE form) that act everywhere uniformly in a spatial domain.

Secondly, employing a cellular dynamics method that includes volume exclusion of cells introduces the possibility that no available sites exist on an inflow boundary when attempting to perform seeding. We refer to such cases as overcrowding. While overcrowding may elucidate problematic ODE model forms and/or parameters, we do not necessarily discard an ODE model or parameter set when overcrowding occurs (though it could inform further ODE model development), neither do we necessarily ignore a seeding event. Rather, we generally describe the inflow of cells as consisting of two stages. In the event of seeding a cell of a type ζ, we first add a cell to a population of nearby cells of type ζ. Nearby cells are not spatially represented, but instead act everywhere homogeneously (i.e., as if they were global) in the spatial domain according to ODE model forms. Before performing all spatial interactions of a simulation step, we attempt to seed each nearby cell. Each nearby cell that is successfully seeded is then removed from the population of nearby cells, and the remaining population of nearby cells act homogeneously in the spatial domain for the simulation step.

Implementation details

The remaining aspects of cellularizing an ODE model depend on the choice in cellular dynamics method to explicitly describe the spatial dynamics that are implicitly represented in an ODE model. Such a choice largely depends on the scale and resolution of the spatial domain, whether on the order of microns, centimeters, or otherwise, which dictates the appropriateness of a particular cellular dynamics method. Various cellular dynamics methods have their own mathematical and computational details that affect the behavior of model cells and their emergent behaviors and properties [35].

To demonstrate cellularization of ODE models in this work, we employed the cellular Potts model (CPM, or Glazier-Graner-Hogeweg model) to model cellular dynamics in analogous spatial models. The CPM represents individual cells and a general medium as deformable, volume-excluding spatial objects in a lattice that defines a spatial domain [27], and is implemented in multicellular modeling and simulation software like CompuCell3D [36], Morpheus [37], and Chaste [38]. Cell motility in the CPM is modeled as the stochastic exchanging of lattice sites at intercellular and cell-medium interfaces. In the CPM, a pair of neighboring lattice sites \( x\in \mathcal{U} \) and \( {x}^{\prime}\in \mathcal{U} \) are randomly selected and the case is considered where the identification at x attempts to copy itself to x, called a copy attempt and denoted σ(x, t) → σ(x, t). A copy attempt occurs with a probability according to a system effective energy \( \mathcal{H} \) that models various cellular properties and behaviors (e.g., volume constraints, adhesion, chemotaxis). For each copy attempt, the change in system effective energy \( \Delta \mathcal{H} \) due to σ(x, t) → σ(x, t) is calculated. Copy attempts that decrease the system effective energy are always accepted, and copy attempts that increase the system effective energy occur according to a Boltzmann acceptance function of the change in system effective energy,

$$ \Pr \left(\sigma \left(x,t\right)\to \sigma \left({x}^{\prime },t\right)\right)={e}^{-\max \left\{0,\frac{\Delta \mathcal{H}}{{\mathcal{H}}^{\ast }}\right\}}. $$
(37)

Here \( {\mathcal{H}}^{\ast } \) is the intrinsic random motility that affects the stochasticity of copy attempts.

In this work, we employed a system effective energy modeling a volume constraint for each cell, adhesion at all interfaces, and chemotaxis,

$$ \mathcal{H}(t)=\sum_s{\lambda}_v\left(\tau \left(s,t\right)\right){\left(v\left(s,t\right)-{v}_c\left(\tau \left(s,t\right)\right)\right)}^2+\sum_{x\in \mathcal{U}}\sum_{x^{\prime}\in \mathcal{N}(x)}\left(1-{\delta}_{\sigma \left(x,t\right),\sigma \left({x}^{\prime },t\right)}\right)J\left(\tau \left(\sigma \left(x,t\right),t\right),\tau \left(\sigma \left({x}^{\prime },t\right),t\right)\right)+\sum_{x\in \mathcal{U}}\sum_c\frac{\lambda_c\left(\tau \left(\sigma \left(x,t\right),t\right)\right)c\left(x,t\right)}{1+{c}_{CM}\left(\sigma \left(x,t\right),t\right)} $$
(38)

The first summation models a volume constraint for each cell with current volume v(s, t) for cell of identification s and type-dependent volume multiplier λv and volume constraint vc. The second summation models adhesion, where δ is the Kronecker-delta, \( \mathcal{N}(x) \) is the neighborhood of a site x in the lattice, and J maps the types occupying the neighboring site pair (x, x) to a contact coefficient. The third summation models logarithmic chemotaxis due to all simulated fields c with type-dependent chemotaxis multiplier λc and measurement cCM of c at the center of mass of a cell occupying a lattice site.

In “Results,” we demonstrate cellularization of two ODE models of viral infection in epithelial cells using the spatial configuration employed in [17]. The first of the two models represents infection in an epithelial sheet, and the second ODE model adds recruitment of immune cells by inflammatory signaling to the first. As such, epithelial cells are arranged in a uniform, planar configuration and fixed in a regular lattice. For simulations without recruitment of immune cells, the spatial domain is two-dimensional. For simulations with recruitment of immune cells, motile immune cells are placed in a second, two-dimensional layer above the epithelial sheet. Boundary conditions were periodic for boundaries parallel to the plane of the sheet, and Neumann for boundaries perpendicular to the plane of the sheet. The planar dimension and lattice size of all simulations is 400 μm and 200 lattice sizes, respectively (lattice site width of 2 μm/site, Table 3). We applied a uniform volume multiplier of 9 and volume constraint of 25 lattice sites (λv and vc in Eq. 38, respectively) to all cells and placed all epithelial cells in a square shape of 5 × 5 × 1 sites. Simulations were performed for 2 weeks to 1 month of simulation time, depending on the ODE model and parameter values, with a simulation step of 5 min. All ODE models considered a total number of epithelial cells NODE equal to 10 million, from which the global cellularization scaling coefficient η was calculated by the number of epithelial cells in the spatial domain as ηNODE= 1600 cells, and the local cellularization scaling coefficient θ was calculated as NODEvcθ= 1 by imposing the cell volume defined in the spatial model on the epithelial cells described in the ODE models. To test the effects of stochasticity, 100 simulation replicas were executed for each spatial model and set of initial conditions and parameters. All simulation measurements were made at intervals of 50 min (ten simulation steps). All simulations were performed in CompuCell3D.

Availability of data and materials

All data generated or analyzed during this study are included in this published article. CompuCell3D scripts to run the simulations of this work are included in the supplementary materials folder “Model.” CompuCell3D is freely available at compucell3d.org and supported on Windows, Mac, and Linux. Datasets generated and analyzed using the CompuCell3D scripts during this work are included in the supplementary materials folder “Results.”

Abbreviations

CPM:

Cellular Potts model

NRMSE:

Normalized root mean squared error

ODE:

Ordinary differential equation

References

  1. Qu Z, Garfinkel A, Weiss JN, Nivala M. Multi-scale modeling in biology: how to bridge the gaps between scales? Prog Biophys Mol Biol. 2011;107(1):21–31. https://doi.org/10.1016/j.pbiomolbio.2011.06.004.

    Article  PubMed  PubMed Central  Google Scholar 

  2. Meier-Schellersheim M, Fraser IDC, Klauschen F. Multiscale modeling for biologists. WIREs Syst Biol Med. 2009;1(1):4–14. https://doi.org/10.1002/wsbm.33.

    Article  CAS  Google Scholar 

  3. Huber F, Schnauß J, Rönicke S, Rauch P, Müller K, Fütterer C, et al. Emergent complexity of the cytoskeleton: from single filaments to tissue. Adv Phys. 2013;62(1):1–112. https://doi.org/10.1080/00018732.2013.771509.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  4. Merks RMH, Glazier JA. A cell-centered approach to developmental biology. Physica A. 2005;352(1):113–30. https://doi.org/10.1016/j.physa.2004.12.028.

    Article  CAS  Google Scholar 

  5. Friedl P, Gilmour D. Collective cell migration in morphogenesis, regeneration and cancer. Nat Rev Mol Cell Biol. 2009;10(7):445–57. https://doi.org/10.1038/nrm2720.

    Article  CAS  PubMed  Google Scholar 

  6. Timm A, Yin J. Kinetics of virus production from single cells. Virology. 2012;424(1):11–7. https://doi.org/10.1016/j.virol.2011.12.005.

    Article  CAS  PubMed  Google Scholar 

  7. Lawson DA, Witte ON. Stem cells in prostate cancer initiation and progression. J Clin Invest. 2007;117(8):2044–50. https://doi.org/10.1172/JCI32810.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Ellis RJ. Macromolecular crowding: obvious but underappreciated. Trends Biochem Sci. 2001;26(10):597–604. https://doi.org/10.1016/S0968-0004(01)01938-7.

    Article  CAS  PubMed  Google Scholar 

  9. Höfling F, Franosch T. Anomalous transport in the crowded world of biological cells. Rep Prog Phys. 2013;76(4):046602. https://doi.org/10.1088/0034-4885/76/4/046602.

    Article  CAS  PubMed  Google Scholar 

  10. Weilandt DR, Hatzimanikatis V. Particle-based simulation reveals macromolecular crowding effects on the Michaelis-Menten mechanism. Biophys J. 2019;117(2):355–68. https://doi.org/10.1016/j.bpj.2019.06.017.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Smith CA, Yates CA. Spatially extended hybrid methods: a review. J R Soc Interface. 2018;15(139):20170931. https://doi.org/10.1098/rsif.2017.0931.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Merks RMH, Perryn ED, Shirinifard A, Glazier JA. Contact-inhibited chemotaxis in de novo and sprouting blood-vessel growth. PLoS Comput Biol. 2008;4(9):e1000163. https://doi.org/10.1371/journal.pcbi.1000163.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Belmonte JM, Clendenon SG, Oliveira GM, Swat MH, Greene EV, Jeyaraman S, et al. Virtual-tissue computer simulations define the roles of cell adhesion and proliferation in the onset of kidney cystic disease. MBoC. 2016;27(22):3673–85. https://doi.org/10.1091/mbc.e16-01-0059.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Sego TJ, Kasacheuski U, Hauersperger D, Tovar A, Moldovan NI. A heuristic computational model of basic cellular processes and oxygenation during spheroid-dependent biofabrication. Biofabrication. 2017;9(2):024104. https://doi.org/10.1088/1758-5090/aa6ed4.

    Article  CAS  PubMed  Google Scholar 

  15. Beauchemin C. Probing the effects of the well-mixed assumption on viral infection dynamics. J Theor Biol. 2006;242(2):464–77. https://doi.org/10.1016/j.jtbi.2006.03.014.

    Article  PubMed  Google Scholar 

  16. Sego TJ, Glazier JA, Tovar A. Unification of aggregate growth models by emergence from cellular and intracellular mechanisms. R Soc Open Sci. 2020;7(8):192148. https://doi.org/10.1098/rsos.192148.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Sego TJ, Aponte-Serrano JO, Gianlupi JF, Heaps SR, Breithaupt K, Brusch L, et al. A modular framework for multiscale, multicellular, spatiotemporal modeling of acute primary viral infection and immune response in epithelial tissues and its application to drug therapy timing and effectiveness. PLoS Comput Biol. 2020;16(12):e1008451. https://doi.org/10.1371/journal.pcbi.1008451.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Murray JM, Goyal A. In silico single cell dynamics of hepatitis B virus infection and clearance. J Theor Biol. 2015;366:91–102. https://doi.org/10.1016/j.jtbi.2014.11.020.

    Article  PubMed  Google Scholar 

  19. Figueredo GP, Siebers P-O, Owen MR, Reps J, Aickelin U. Comparing stochastic differential equations and agent-based modelling and simulation for early-stage cancer. PLoS ONE. 2014;9(4):e95150. https://doi.org/10.1371/journal.pone.0095150.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Glont M, Arankalle C, Tiwari K, Nguyen TVN, Hermjakob H, Malik-Sheriff RS. BioModels Parameters: a treasure trove of parameter values from published systems biology models. Bioinformatics. 2020;36(17):4649–54. https://doi.org/10.1093/bioinformatics/btaa560.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Glont M, Nguyen TVN, Graesslin M, Hälke R, Ali R, Schramm J, et al. BioModels: expanding horizons to include more modelling approaches and formats. Nucleic Acids Res. 2018;46(D1):D1248–53. https://doi.org/10.1093/nar/gkx1023.

    Article  CAS  PubMed  Google Scholar 

  22. Gallagher ME, Brooke CB, Ke R, Koelle K. Causes and consequences of spatial within-host viral spread. Viruses. 2018;10(11):627. https://doi.org/10.3390/v10110627.

    Article  CAS  PubMed Central  Google Scholar 

  23. Michael Lavigne G, Russell H, Sherry B, Ke R. Autocrine and paracrine interferon signalling as ‘ring vaccination’ and ‘contact tracing’ strategies to suppress virus infection in a host. Proc R Soc B Biol Sci. 2021;288(1945):20203002. https://doi.org/10.1098/rspb.2020.3002.

    Article  CAS  Google Scholar 

  24. Sun J, Vera JC, Drnevich J, Lin YT, Ke R, Brooke CB. Single cell heterogeneity in influenza A virus gene expression shapes the innate antiviral response to infection. PLoS Pathog. 2020;16(7):e1008671. https://doi.org/10.1371/journal.ppat.1008671.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Price I, Mochan-Keef ED, Swigon D, Ermentrout GB, Lukens S, Toapanta FR, et al. The inflammatory response to influenza A virus (H1N1): an experimental and mathematical study. J Theor Biol. 2015;374:83–93. https://doi.org/10.1016/j.jtbi.2015.03.017.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Best K, Perelson AS. Mathematical modeling of within host Zika virus dynamics. Immunol Rev. 2018;285(1):81–96. https://doi.org/10.1111/imr.12687.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Graner F, Glazier JA. Simulation of biological cell sorting using a two-dimensional extended Potts model. Phys Rev Lett. 1992;69(13):2013–6. https://doi.org/10.1103/PhysRevLett.69.2013.

    Article  CAS  PubMed  Google Scholar 

  28. Ball SC, Abraha A, Collins KR, Marozsan AJ, Baird H, Quiñones-Mateu ME, et al. Comparing the ex vivo fitness of CCR5-Tropic human immunodeficiency virus type 1 isolates of subtypes B and C. JVI. 2003;77(2):1021–38. https://doi.org/10.1128/JVI.77.2.1021-1038.2003.

    Article  CAS  Google Scholar 

  29. Quiñones-Mateu ME, Ball SC, Marozsan AJ, Torre VS, Albright JL, Vanham G, et al. A dual infection/competition assay shows a correlation between ex vivo human immunodeficiency virus type 1 fitness and disease progression. J Virol. 2000;74(19):9222–33. https://doi.org/10.1128/JVI.74.19.9222-9233.2000.

    Article  PubMed  PubMed Central  Google Scholar 

  30. Olmsted SS, Padgett JL, Yudin AI, Whaley KJ, Moench TR, Cone RA. Diffusion of macromolecules and virus-like particles in human cervical mucus. Biophys J. 2001;81(4):1930–7. https://doi.org/10.1016/S0006-3495(01)75844-4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Leal J, Smyth HDC, Ghosh D. Physicochemical properties of mucus and their impact on transmucosal drug delivery. Int J Pharm. 2017;532(1):555–72. https://doi.org/10.1016/j.ijpharm.2017.09.018.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Krummel MF, Bartumeus F, Gérard A. T cell migration, search strategies and mechanisms. Nat Rev Immunol. 2016;16(3):193–201. https://doi.org/10.1038/nri.2015.16.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Pearson JE, Krapivsky P, Perelson AS. Stochastic theory of early viral infection: continuous versus burst production of virions. PLoS Comput Biol. 2011;7(2):e1001058. https://doi.org/10.1371/journal.pcbi.1001058.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Sluka JP, Fu X, Swat M, Belmonte JM, Cosmanescu A, Clendenon SG, et al. A liver-centric multiscale modeling framework for xenobiotics. PLoS ONE. 2016;11(9):e0162428. https://doi.org/10.1371/journal.pone.0162428.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Osborne JM, Fletcher AG, Pitt-Francis JM, Maini PK, Gavaghan DJ. Comparing individual-based approaches to modelling the self-organization of multicellular tissues. PLoS Comput Biol. 2017;13(2):e1005387. https://doi.org/10.1371/journal.pcbi.1005387.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Swat MH, Thomas GL, Belmonte JM, Shirinifard A, Hmeljak D, Glazier JA. Multi-scale modeling of tissues using CompuCell3D. Methods Cell Biol. 2012;110:325–66. https://doi.org/10.1016/B978-0-12-388403-9.00013-8.

    Article  PubMed  PubMed Central  Google Scholar 

  37. Starruß J, de Back W, Brusch L, Deutsch A. Morpheus: a user-friendly modeling environment for multiscale and multicellular systems biology. Bioinformatics. 2014;30(9):1331–2. https://doi.org/10.1093/bioinformatics/btt772.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Cooper FR, Baker RE, Bernabeu MO, Bordas R, Bowler L, Bueno-Orovio A, et al. Chaste: cancer, heart and soft tissue environment. J Open Source Softw. 2020;5(47):1848. https://doi.org/10.21105/joss.01848.

    Article  Google Scholar 

Download references

Acknowledgements

This research was supported in part by Lilly Endowment, Inc., through its support for the Indiana University Pervasive Technology Institute. The funders had no role in manuscript preparation or the decision to submit the work for publication.

Funding

The authors acknowledge funding from National Institutes of Health grants U24 EB028887 and R01 GM122424 and National Science Foundation grant NSF 1720625.

Author information

Authors and Affiliations

Authors

Contributions

TJS, JAS, JFG, and JAG contributed to conceptualization and methodology, and review and editing of the manuscript, of this work. TJS and JAG administered and provided resources to accomplish this work. TJS performed the investigation, wrote all software, performed validation, formal analysis, data curation, and visualization, and wrote the initial draft of the manuscript, of this work. JAG supervised and acquired the funding to accomplish this work. All authors read and approved the final manuscript.

Corresponding author

Correspondence to T. J. Sego.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

JAG is the owner/operator of Virtual Tissues for Health, LLC, which develops applications of multiscale tissue models in medical applications.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1.

CompuCell3D scripts and Datasets generated and analyzed.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Sego, T.J., Aponte-Serrano, J.O., Gianlupi, J.F. et al. Generation of multicellular spatiotemporal models of population dynamics from ordinary differential equations, with applications in viral infection. BMC Biol 19, 196 (2021). https://doi.org/10.1186/s12915-021-01115-z

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12915-021-01115-z

Keywords