Multiscale mathematical models for biological systems

Abstract

Life activities are extremely complex phenomena in nature. From molecular signaling regulation to multi-cellular tissue formation and so on, the biological system consists of multiple temporal, spatial and functional scales. Multiscale mathematical models have extensive applications in life science research due to their capacity of appropriately simulating the complex multiscale biological systems. Many mathematical methods, including deterministic methods, stochastic methods as well as discrete or rule-based methods, have been widely used for modeling biological systems. However, the models at single scale are not sufficient to simulate complex biological systems. Therefore, in this paper we give a survey of two multiscale modeling approaches for biological systems. One approach is continuous stochastic method that couples ordinary differential equations and stochastic differential equations; Another approach is hybrid continuous-discrete method that couples agent-based model with partial differential equations. We then introduce the applications of these multiscale modeling approaches in systems biology and look ahead to the future research.

Key words： Multi-scale model; partial differential equations; agent-based model; systems biology

1 Introduction

Life activity is one of the most complex natural phenomena. Biological systems are constituted of multi-functional networks that regulate different temporal states and spatial domains of organism, maintaining growth, development, and regeneration of the organism. Ranging from elementary amino acid molecules determining protein functions to cell-cell interactions regulating secretion of hormones, the biological systems are multiscale in time, space and function [13]. Mathematical models that quantitatively describe and integrate different biological scales can help understand, predict and control collective behavior of biological systems and provide guidance for further experiments. Therefore, multiscale models are uniquely adequate and powerful tools for biological systems modeling.

Multiscale modeling plays important roles and has important applications in life science. The 2013 Nobel Prize in Chemistry was awarded to Martin Karplus, Michael Levitt and Arieh Warshel, three computational scientists working in the life sciences, for their contributions to the development of multiscale models of complex chemical systems. They built theoretical models based on quantum mechanics, classical mechanics and hybrid quantum-classical mechanics to study the motion of protein molecules and enzyme-catalyzed reaction mechanism by means of computational simulation. The tradition of modern science is the combination of positivism and rationalism. Positivism refers to experiments, data, observations, etc. Rationalism refers to logic, deduction, especially mathematics. In fact, the combination of mathematical models and experimental studies was vividly reflected in the development of physics in the 20th century, so that physics has achieved rich results and considerable development, among which relativity and quantum mechanics are typical cases of this combination. Similarly, modern biology has gradually developed from the traditional experimental approach to its combination with quantitative analysis, and from the macroscopic morphology and phenotype of organs and tissues to the in-depth study of molecular and cellular mechanisms. In the past decades, a wealth of experimental observation and knowledge has accumulated. Meanwhile, more and more technologies to generate high-throughput data, such as a variety of omics data, have been developed. Utilizing mathematical and computational tools to mine and integrate these data has become an effective way for biological research.

The behavior, functions and states of biological systems are determined by the collective interactions of their constituent parts. As such, it’s more reasonable to study the biological system as a whole, rather than focusing only on individual genes, proteins, or cells, for better understanding of regulatory mechanisms and evolutionary processes. In view of this, systems biology tries to integrate different levels of information to understand how biological systems perform their functions. It is the discipline that studies the architecture of the biological system, interactions between the internal components and evolutionary dynamics of the system under various internal and external conditions. It is a frontier field of life science, and also an interdiscipline of biology, mathematics and systems science [2, 28].

Most biological systems are dynamic and stochastic, involving different temporal, spatial and functional scales. For instance, biochemical reactions associated with protein-protein interactions belong to short-time scales (seconds/minute), while the developmental changes of phenotypes such as cell differentiation and apoptosis are typically on a long-time scale (hour/day). In terms of spatial scale, from protein size (nanometers) to cell diameter (microns) to tissue size (millimeters) and to human body (meters), it shows a step increase from micro scale to macro scale [7]. Obviously, it is difficult to fully investigate multi-scale biological systems through experimental methods alone. Alternatively, it is necessary to develop mathematical models and computational simulations to more comprehensively study the integrative behavior and dynamic properties of biological system. Therefore, multi-scale mathematical model plays an indispensable role in system biology modeling.

The multi-scale model is used to simulate a biological system with multiple interacting levels (e.g., molecular level, subcellular, cellular level, microenvironmental level, tissue level, etc.) that the single-scale model can not be competent. Moreover, the multi-scale model enables information exchange between multiple levels so it can provide more information for each level than the single-scale model. To accurately describe complex, multiscale biological systems and make appropriate predictions, a good mathematical model should consider the main factors affecting the system as comprehensively as possible, but at the same time, it should reduce complexity to simplicity as much as possible. Just as Albert Einstein said, “Make everything as simple as possible, but no simpler.” In this sense, mathematical modeling is not only a kind of technology, but also an art for scientific interpretation and discovery with mathematical language.

2 Common mathematical methods for modeling biological systems

Many mathematical methods have been used for simulating biological systems, including deterministic methods, stochastic methods, and rule-based methods, etc. Here we introduce those different approaches and their applications in systems biology.

2.1 Differential equations

The deterministic dynamical systems are often described by using differential equations, which has a large number of instances and wide applications in biological modeling. A typical example is molecular dynamics simulation [7, 25] that simulate the motion of atoms or molecules according to Newton’s second law of motion as follows,

(1)m_i\frac{\mathrm{d}{v}_i}{\mathrm{d}t}=F_i,

where m_i is the mass of the i-th molecule, v_i is its velocity and F_i is the total external force acting on the molecule. Hence, the motion of molecules can be simulated provided F_i is determined. Although this method has been widely applied in many biological questions, yet it’s hard to use this method to simulate systems beyond microsecond scale due to the heavy computational requirement [30].

Many biological processes, such as biological oscillations, often occur at multiple time scales much larger than microsecond scale (e.g., 1-second heart beats and 24-hour cell cycles). For such systems, it’s impossible to simulate molecular dynamics at atomic level. Instead, ordinary differential equations (ODEs) are often employed to simulate the macroscopic properties of these systems. For example, based on the law of mass action [41], ODEs model could be constructed to simulate biochemical reactions (e.g., gene regulation, protein-protein interaction, metabolic reaction, etc.). For a simple reversible reaction, phosphorylation and dephosphorylation between protein A and B, its chemical reaction formula is

\mathrm{A}\rightleftarrows \beta \alpha \mathrm{B},

where α and β represent the phosphorylation rate and dephosphorylation rate of A and B, respectively. If concentrations of A and B are denoted as x and y, respectively, then the production rate of B is \alpha x, and its dephosphorylation rate is \beta y. Hence, the change rate of the concentration of B is

(2)\frac{\mathrm{d}y}{\mathrm{d}t}=\alpha x- \beta y.

The method of deriving ODEs based on the law of mass action is often applied to modeling the dynamics of signaling transduction or pathways [1].

If we consider the diffusion of molecules in space, then we can use partial differential equation (PDE), e.g., reaction-diffusion equation, to describe spatio-temporal distribution and evolution of molecule concentration, which can be derived according to the law of conservation of mass. For example, in the above case, if protein B diffuses in the space, then according to Fick's first law, the change of its concentration y within the unit time can be described as follows,

(3)J=-D\mathrm{\nabla }y,

where D is the diffusion coefficient. The above equation reflects that the molecules tend to move along the direction from high concentration to low concentration. On the other hand, the changing rate of protein B owing to phosphorylation and dephosphorylation is R=\alpha x- \beta y, as mentioned above. According to the law of conservation of mass or Fick's second law, the total changing rate of concentration y is

(4)\frac{\mathrm{\partial }y}{\mathrm{\partial }t}=-\mathrm{d}\mathrm{i}\mathrm{v}(J)+ R,

where div represents divergence operator, \mathrm{d}\mathrm{i}\mathrm{v}\left(J\right) represents the density of flux of vector field J at each location. As such, \mathrm{d}\mathrm{i}\mathrm{v}\left(J\right)>0 means that the given location is the source of J, and indicates that there is a sink of J at the given location. If the diffusion coefficient D is isotropous, then we get the following reaction-diffusion equation that describes the spatio-temporal evolution of the protein concentration y,

(5)\frac{\mathrm{\partial }y}{\mathrm{\partial }t}=D\mathrm{\Delta } y+\alpha x-\beta y.

The reaction-diffusion equation described above is a powerful and concise mathematical tool to analyze the dynamical mechanisms of pattern formation [15, 19], an interesting phenomenon in biology and ecology, which was firstly studied by the British mathematician Alan Mathison Turing.

2.2 Stochastic methods

Randomness plays an important role in signaling regulation and cell functioning due to thermodynamic fluctuations of molecules’ accumulation or local migration/immigration of cells. As such, stochastic methods are increasingly gaining traction in biological systems modeling.

Markov chain can be used to represent biochemical processes. Based on the probabilistic meaning of transition probability matrix of Markov chain, Gillespie algorithm [20, 21] proposed by Daniel Thomas Gillespie is an effective stochastic simulation algorithm for biochemical reactions.

Chemical master equation is another method for modeling stochastic chemical reaction system [37, 45]. For example, for a reaction involving N+1 molecules, there are 0,1,2 ,\dots ,N molecules of protein B and accordingly N,N-1,N-2 ,\dots ,0 molecules of protein A. Assume that the probability of having n molecules of B at time t is P(n, t), then the changing rate of this probability is constituted of four transitions of corresponding states. Therefore, we have N+1 ODEs for N+1 molecules involved in the reaction,

(6)\frac{\mathrm{d}P(n,t)}{\mathrm{d}t}=\alpha (N-n+1)P( n-1,t)+\beta (n+1 )P(n+1,t)-\beta nP( n,t) -\alpha ( N-n)P(n ,t).

Given the initial distribution P(n, 0)(n=0,1,\dots ,N), the above ODEs can be solved to obtain the time-dependent probability distribution P(n,t).

When the number of molecules is relatively large, the master equation can be approximated with Langevin equation, a type of stochastic differential equation (SDE). Langevin equations can be used to simulate the stochastic dynamics of the signaling pathways with extrinsic or intrinsic noises [64]. It has the following form:

(7)\mathrm{d}X=f(X ,t)\mathrm{d}t+g(X ,t)\mathrm{d}W,

where X is a stochastic variable representing the number of particles. The first term in the above equation is the drift term, and the second term is the diffusion term. W is Wiener process or also called as standard Brown motion.

The probability distribution of the stochastic system can also be described using Fokker−Plank equation [29], a type of PDE. The Fokker−Plank equation corresponding to the above Eq. (7) is

(8)\frac{\mathrm{\partial }P(x,t)}{\mathrm{\partial }t}=\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{x}^2}{g}^2(x, t)P(x,t )-\frac{\mathrm{\partial }}{\mathrm{\partial }x} \left[f(x,t)+g^{\prime }(x,t)g(x,t)\right]P(x ,t).

Fokker−Plank equation can be used to conveniently analyze the time evolution and steady-state of the probability distribution. However, the analytical derivation and solution of Fokker−Plank equation will become difficult when the system is complex or the original equations are more coupled. In this situation, Monte Carlo method is often used to simulate the probability distribution of the system.

2.3 Discrete or rule-based methods

Although continuous methods such as differential equations can succinctly describe the dynamic changes in the quantity of some objects in biological systems, due to the complexity of biological systems, the continuity hypothesis may no longer applicable to some processes or systems, so it may not be sufficiently adequate to use continuous functions to directly describe those properties. Therefore, some discrete modeling methods must be used for the modeling of many biological systems, such as Boolean network model [27], difference equation model, Cellular automata model [11] or agent-based model (ABM) [6].

3 Two types of multiscale mathematical models

The above various types of methods are mainly used to describe one certain scale/level of the biological system. As mentioned above, single-level model is not sufficient to characterize complex multi-scale biological systems. In view of this, we here introduce two kinds of multiscale modeling approaches. One is continuous method, i.e., multiscale modeling by coupling ODEs and SDEs, called as multiscale SDE model; another is continuous-discrete hybrid method, i.e., multiscale systems modeling by coupling PDEs and ABM, called as PDE-Agent model.

3.1 Multiscale SDE model

A large number of biological systems are dynamic and stochastic. How to consider effects of random factors across different scales on the dynamics of the system is an interesting problem. The multiscale SDE model can be used to simulate the dynamic properties of such stochastic systems, for instance, the dynamic changes of signaling pathways and cell population dynamics. It employs ODEs to describe biochemical reactions such as phosphorylation and dephosphorylation involved in the signaling pathways and SDEs to simulate dynamic changes in cell number. In this way, it takes into account random factors such as stochastic switching of cellular phenotypes regulated by signaling pathways or local diffusion of the cell population.

3.1.1 Molecular scale: signaling pathway

As mentioned above, the time-change of concentrations of components (e.g., proteins) in the signaling pathways are often described using ODEs. As described by Eq. (2), the ODEs model can be built according to the law of Mass action. However, the detailed mechanisms of biochemical reactions involved in the signaling pathways are often not fully clear, or the sample of experimental data is not large enough to estimate the complete set of parameters in the mass-action kinetics. Alternatively, Hill functions [12, 31] are often used to simulate the dynamics of signaling pathways based on Michealis-Menten kinetics. For example, if protein A can activate protein B, then we can use the following equation to describe this regulation:

(9)\frac{\mathrm{d}y}{\mathrm{d}t}=\frac{V\cdot x^n}{K+x^n}-d\cdot y,

where x and y represent concentrations of proteins A and B, respectively. V is the maximal activation rate of B promoted by A, K is the Michealis-Menten constant, and n is Hill coefficient. d is the degradation rate of activated protein B.

3.1.2 Cellular scale: dynamics of cell populations

We here use a set of SDEs [46] to simulate dynamic changes in cell numbers of different interacting cell populations. For example, the changes in numbers of pre-osteoblast (OBp) and active osteoblast (OBa) during bone remodeling can be described using the following equations:

(10)\mathrm{d}{X}_{\mathrm{O}\mathrm{B}\mathrm{p}}={\tilde{D}}_{\mathrm{M}\mathrm{S}\mathrm{C}}\cdot n_{\mathrm{M}\mathrm{S}\mathrm{C}}\mathrm{d}t+(P_{\mathrm{O}\mathrm{B}\mathrm{p}}-{\tilde{D}}_{\mathrm{O}\mathrm{B}\mathrm{p}}-d_{\mathrm{O}\mathrm{B}\mathrm{p}}) \cdot X_{\mathrm{O}\mathrm{B}\mathrm{p}}\mathrm{d}t+X_{\mathrm{O}\mathrm{B}\mathrm{p}}\cdot {\sigma }_1\mathrm{d}{W}_1,

(11)\mathrm{d}{X}_{\mathrm{O}\mathrm{B}\mathrm{a}}={\tilde{D}}_{\mathrm{O}\mathrm{B}\mathrm{p}}\cdot X_{\mathrm{O}\mathrm{B}\mathrm{p}}\mathrm{d}t+(P_{\mathrm{O}\mathrm{B}\mathrm{a}}-d_{\mathrm{O}\mathrm{B}\mathrm{a}}) \cdot X_{\mathrm{O}\mathrm{B}\mathrm{a}}\mathrm{d}t+X_{\mathrm{O}\mathrm{B}\mathrm{a}}\cdot {\sigma }_2\mathrm{d}{W}_2,

where X_{\mathrm{O}\mathrm{B}\mathrm{p}} and X_{\mathrm{O}\mathrm{B}\mathrm{a}} are stochastic processes representing cell densities of OBp and OBa, respectively. {\tilde{D}}_{\mathrm{M}\mathrm{S}\mathrm{C}} and {\tilde{D}}_{\mathrm{O}\mathrm{B}\mathrm{p}} respectively represent differentiation rate of Mesenchymal stem cell (MSC) to OBp and differentiation rate of OBp to OBa. P_{\mathrm{O}\mathrm{B}\mathrm{p}} and P_{\mathrm{O}\mathrm{B}\mathrm{a}} represent self-proliferation rates of OBp and OBa, respectively. d_{\mathrm{O}\mathrm{B}\mathrm{p}} and d_{\mathrm{O}\mathrm{B}\mathrm{a}} are apoptosis rates of the respective cell types. W_i(i=1,2) is Wiener process or also called as standard Brown motion, which is used here to describe the local migration/immigration and stochastic regulation of cell activities by signaling pathways. Accordingly, {\sigma }_i(i=1 ,2) represents the diffusion coefficient. The above SDEs can be numerically solved using Euler-Maruyama method [23].

In addition, it can be assumed that the binding and disassociation of proteins in the signaling pathways are rather faster than the switching of cell phenotypes and changes in densities of cell populations [34, 48]. More specifically, the signaling pathway is at shorter time scale (minutes/hours), while the cell activity is at longer time scale (days). Therefore, the quasi-steady-state of ODEs for the signaling pathways can be derived, which can be further connected to cellular scale to regulate numbers of different types of cells via, for example, {\tilde{D}}_{\mathrm{M}\mathrm{S}\mathrm{C}} and {\tilde{D}}_{\mathrm{O}\mathrm{B}\mathrm{p}}. In this way, the two biological levels with different time scales can be linked together.

Note that this type of model does not explicitly represent the stochastic factors in the signaling pathways. Instead, the stochastic factors are implicated in the SDEs describing cellular scale. In this way, it significantly reduces the model complexity, the number of unknown parameters and computation heavy as well, without loss of consideration of the stochastic factors in the cell activities regulated by signaling pathways.

3.2 Multiscale PDE-Agent model

Although the continuous methods can simulate dynamics of signaling pathways and cell populations, but this approach falls short on simulating interactions between individual cells and cell-microenvironment interaction. To more precisely model spatio-temporal evolution of some complex biological systems, we here introduce a kind of hybrid continuous-discrete model that combines a continuous module that uses PDEs to describe microenvironment and a discrete module that employs ABM to simulate cell activity.

Agent-based modeling (ABM) is a type of rule-based modeling strategy that borrows ideas of Agent, a concept in artificial intelligence. ABM has wide applications in many complex biological systems and social networks. In multi-cellular systems, each cell can be viewed as an intelligent agent that receives stimuli from the microenvironment and responds to them. Each agent is assigned with a set of extrinsic and intrinsic rules to evolve. These rules always have random characteristics. Hence, a large amount of single cell activities constitutes the collective behavior of cell population, often leading to emergency properties.

The rules in ABM are always defined in a two- or three-dimensional lattice. The whole lattice represents a slide or a portion of a tissue. The grid interval is set to diameter size of typical cells. Each cell occupies a grid in the lattice.

3.2.1 Molecular scale: signaling pathway

As described above, the intracellular signaling pathways are often modeled by a set of ODEs. We omit the details here.

3.2.2 Cellular scale: cell activity

Cell activities (e.g., migration, differentiation, proliferation and apoptosis) can be simulated using a set of rules. These rules depend on some key signaling proteins or transcriptional factors (TFs) downstream of the signaling pathway. For example, TFs Runx2 and Osx can promote differentiation of MSC to OBp [22, 33, 35]. So the probability of differentiation of MSC to OBp ( P_{\mathrm{M}\mathrm{S}\mathrm{C}\to \mathrm{O}\mathrm{B}\mathrm{p}}^{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}) depends on levels of activated Runx2 and Osx, which can be modeled as follows,

(12)P_{\mathrm{M}\mathrm{S}\mathrm{C}\to \mathrm{O}\mathrm{B}\mathrm{p}}^{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}=\left(\frac{V_{D1,\mathrm{R}\mathrm{u}\mathrm{n}\mathrm{x}2}\cdot [\mathrm{R}\mathrm{u}\mathrm{n}\mathrm{x}2]}{K_{D1,\mathrm{R}\mathrm{u}\mathrm{n}\mathrm{x}2}+[\mathrm{R}\mathrm{u}\mathrm{n}\mathrm{x}2]}+\frac{V_{D1,\mathrm{O}\mathrm{s}\mathrm{x}}\cdot [\mathrm{O}\mathrm{s}\mathrm{x}]}{K_{D1,\mathrm{O}\mathrm{s}\mathrm{x}}+[\mathrm{O}\mathrm{s}\mathrm{x}]}\right)\cdot p_{\mathrm{M}\mathrm{S}\mathrm{C}\to \mathrm{O}\mathrm{B}\mathrm{p}}^{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}},

where p_{\mathrm{M}\mathrm{S}\mathrm{C}\to \mathrm{O}\mathrm{B}\mathrm{p}}^{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}} represents basal rate of MSC differentiating to OBp. V_{D1,\mathrm{R}\mathrm{u}\mathrm{n}\mathrm{x}2} and V_{D1,\mathrm{O}\mathrm{s}\mathrm{x}} are Hill regulatory factors of Runx2 and Osx, respectively. K_{D1,\mathrm{R}\mathrm{u}\mathrm{n}\mathrm{x}2} and K_{D1,\mathrm{O}\mathrm{s}\mathrm{x}} are Hill regulatory coefficients of Runx2 and Osx, respectively [61].

3.2.3 Microenvironment scale: extracellular environment

Various types of biochemical molecules (e.g., glucose, growth factors, and cytokines) diffuse within the extracellular microenvironment. The secretion, diffusion, permeation and uptake of those molecules can be described by a set of reaction-diffusion equations. For example, tumor cells secrete vascular endothelial growth factor (VEGF) to induce angiogenesis. After secretion, VEGF diffuses in the microenvironment and can be uptake by endothelial cells [56]. The above processes can be described by the following reaction-diffusion equations:

(13)\frac{\mathrm{\partial }V}{\mathrm{\partial }t}=D_V\mathrm{\Delta }V+{\chi }_{\mathrm{t}\mathrm{u}\mathrm{m}}(t,x)S_V-{\chi }_{\mathrm{v}\mathrm{e}\mathrm{s}}(t,x)q_{V}V-{\delta }_{V}V,

where V is VEGF concentration, and D_V is its diffusion coefficient. S_V is the secretion rate of VEGF by tumor cells. q_V is the uptake rate of VEGF by endothelial cells. {\delta }_V is the degradation rate of VEGF. {\chi }_{\mathrm{t}\mathrm{u}\mathrm{m}}(t, x) represents time-dependent characteristic function, which equals to 1 in the tumor region and 0 in other locations. Similarly, {\chi }_{\mathrm{v}\mathrm{e}\mathrm{s}}(t,x) equals to 1 when there is a vessel at x and 0 otherwise. During each simulation step, χ_tum and χ_ves are updated according to the spatial distributions of tumor cells and blood vessels.

Given appropriate initial condition and boundary condition, the above reaction-diffusion equations can be numerically solved using, for instance, finite difference method [42].

3.2.4 Tissue scale: tissue growth and angiogenesis

The tissue scale directly connects exterior environment to the interior cells. The exterior stimuli (such as mechanical stimulation) are delivered to cells through the tissue. On the other hand, functional cells promote or influence the tissue growth. For example, in the bone remodeling, Osteoblasts and osteoclasts promote bone formation and absorption, respectively. The balance between these two processes controls the process of bone remodeling, while the external mechanical loading influences the accumulation of bone mass. Therefore, tissue growth depends on the functional activities of various types of cells and the stimulation of the external environment.

Angiogenesis is a critical process in many pathological or physiological phenomena. Thus, it is sometimes one of the main objects of concern when developing models at tissue scale. It is assumed that endothelial cells located at the tip of capillary buds control the movement of the blood vessel, and that the chemotaxis of these cells with respect to VEGF and the haptotaxis to fibronectin are the main factors affecting their movement [56]. Therefore, the probability of migration of tip endothelial cells can be defined as follows (Fig.1),

Fig.1 Illustration of rules of migration of tip endothelial cells.

Full size|PPT slide

(14)P_k\propto \left(\frac{\alpha k_V}{k_V+V}\mathrm{\nabla }V+\lambda \mathrm{\nabla }F\right)\cdot l_k,k=1,2,3,4,

where V is VEGF concentration, and F is fibronectin concentration. l_k is a unit vector along the k-th direction. \frac{\alpha k_V}{k_V+V}\mathrm{\nabla }V describes the chemotaxis of the tip endothelial cells to the gradient of VEGF concentration [4]. α is the chemotaxis coefficient. k_V is a constant controlling the weight of VEGF concentration in the chemotaxis sensitivity. \lambda \mathrm{\nabla }F represents the extent of haptotaxis of endothelial cells to fibronectin, with λ as haptotaxis coefficient.

The probability of tip endothelial cell remaining stationary can be defined as the average of P₁, P₂, P₃ and P₄. After normalization, these probability gives possibilities of the tip endothelial cell moving up, down, right, left or keeping stationary.

The algorithm for simulating angiogenesis is as follows:

(1) Calculate migration probability of the endothelial cell

a. Numerically solve the PDEs of microenvironmental factors, such as VEGF and fibronectin (e.g., Eq. (13)). Compute P₁−P₄ and P₅ according to Equation (14);

b. Normalize the probability: {\tilde{P}}_i=P_i/{\sum }_{j=1}^5{P}_j, i=1,2,\dots ,5; Define interval I_1= (0,{\tilde{P}}_1], I_i=( {\sum }_{j=1}^{i-1}{\tilde{P}}_j, {\sum }_{j=1}^i{\tilde{P}}_j], i=2,\dots ,5.

(2) For each sprout tip cell, check whether the age of the corresponding vessel is greater than 18 and whether enough space in the neighbor grids is available.

a. Vessel branching: if the above conditions are satisfied, generate two random numbers between 0 and 1 and then perform simulation of vessel branching. For example, If r₁ \in I₂, r₂ \in I₃, then move two endothelial cells down and right, respectively; if random numbers are located within other intervals, perform similar procedures.

b. Sprout migration: if the above conditions are not satisfied, generate another random number r between 0 and 1 and then perform simulation of vessel migration. For example, If r \in I₃, then move the tip endothelial cell right; if random numbers are located within other intervals, perform similar procedures.

(3) Vessel fusion: if two vessels meet each other, then they merge into a single blood vessel and continue to grow.

3.2.5 Integration of different scales

The algorithm for connecting different scales during each step in the simulation is summarized as follows (Fig.2).

Fig.2 Flowchart of computational algorithm for multiscale PDE-Agent model.

Full size|PPT slide

(1) At microenvironmental scale, numerically solve PDEs to get the spatio-temporal distributions of concentrations of microenvironmental factors (e.g., growth factor, glucose, oxygen, and VEGF);

(2) At molecular scale, use the above solutions of microenvironmental factors as input for the ODEs of the intracellular signaling pathways;

(3) At cellular scale, simulate cell activities (e.g., migration, differentiation, proliferation, and apoptosis) according to the defined rules that depend on the status of the signaling pathways and surrounding environment, and update χ_tum(t,x) according to the profile of tumor growth;

(4) At tissue scale, simulate migration of tip endothelial cells and branching of vessel buds, and update χ_ves(t,x) according to the profile of the remodeled vasculature network;

(5) Return to the step (1), substitute the updated χ_tum(t,x) and χ_ves(t,x) into PDEs for next round of numerical simulation.

4 Parameter estimation

How to estimate key parameters in the above model is an important question. For example, in the signaling pathway modeling, how to infer interaction network of proteins from biochemical experimental data (e.g., Western blotting) or proteomics data (e.g., protein mass spectrometry data, reverse protein phosphorylation array (RPPA) data) and to estimate parameters in the ODEs model describing dynamics of signaling pathways, is a substantially important step for multiscale modeling. Only by establishing a dynamic network model at the molecular level and estimating its parameters from experimental data can we simulate cell activity and tissue growth on the basis of molecular mechanisms.

The experimental data for the signaling pathways reflect the degree of activities of related proteins (e.g., protein phosphorylation) under different conditions or at different time points under certain experimental conditions, which are actually “results” of the dynamic changes of signaling pathways. The differential equation model and parameter estimation attempt to explore this “process” of dynamic changes from the information provided by these experimental results. Therefore, in a broad sense, the modeling and parameter estimation of signaling network actually require solving an “inverse problem”. However, the inverse problem is a challenging problem that has not been fully addressed so far. The inverse problem often requires concrete analysis of specific problems and it lacks a unified solution method. In practical problems, it is necessary to mine the prior information of the solution as much as possible, and impose some constraints according to the background and requirements of the problem, so as to narrow the search space of the solution and make it easier to obtain the solution in line with the actual meaning of the problem.

Nonlinear least square objective function can be defined as follows for the estimation of unknown parameters in the ODE models,

(15)\hat{\theta }=\underset{\theta \in \mathrm{\Theta }}{\arg min} \sum _{i=1}^m\sum _{l=1}^n{({y_i}^{\mathrm{s}\mathrm{i}\mathrm{m}}(t_l,\theta )-{y_i}^{\exp }(t_l,\theta ))}^2 ,

where {y_i}^{\mathrm{s}\mathrm{i}\mathrm{m}}( t_l,\theta ) and {y_i}^{\exp }(t_l,\theta ) are the simulated and experimental data of concentrations of the i-th protein under the l-th condition (or at the l-th time point t_l) with parameters \theta =({\theta }_1,{\theta }_2,\dots ,{\theta }_N). \mathrm{\Theta } \subset R^N is the parameter space. The searching interval for each parameter is set to a range according to the experimental observations or the biochemical meaning of Michaelis-Menten kinetics. Genetic algorithm [26] can be used for optimization of the objective function in the Eq. (15).

Here we briefing the procedure of using the genetic algorithm to estimate parameters.

(1) Set the searching space of parameters (\mathrm{\Theta } ) according to prior information and randomly select a set of initial parameters.

(2) Repeat the following steps：

a. for each set of parameters, numerically solve ODEs using 4th order Runge-Kutta method;

b. evaluate the fitness of each parameter set based on Eq. (15)−the smaller the objective function value, the greater the fitness score;

c. select excellent individuals according to their fitness scores (those with high scores are more likely to be selected);

d. Change the parameter set via crossover and mutation.

(3) check whether the stopping criterion (e.g., the objective function value no longer decreasing) is satisfied. If yes, then output the solution; otherwise, return to the step (2).

5 Sensitivity analysis and model robustness

The sensitivity analysis is often employed to examine the extent of change of the system in response to the small perturbations of parameters. The relative sensitivity coefficient [43] of variable y_i (i=1,\dots ,m) at time point t\in [0,T ] with respect to {\theta }_j (j=1,\dots ,N) can be calculated via the following formula,

(16)S_{ij}(t)=\frac{\mathrm{\partial }{y}_i/\mathrm{\partial }{\theta }_j}{y_i/{\theta }_j}.

For small \mathrm{\Delta }{\theta }_j (e.g., \mathrm{\Delta }{\theta }_j/\phantom{\mathrm{\Delta }{\theta }_j{\theta }_j}{\theta }_j = 1%),

(17)S_{ij}(t)\approx \frac{\mathrm{\Delta }{y}_i}{y_i}/\phantom{\frac{\mathrm{\Delta }{y}_i}{y_i}\frac{\mathrm{\Delta }{\theta }_j}{{\theta }_j}}\frac{\mathrm{\Delta }{\theta }_j}{{\theta }_j}=\frac{y_i(t,{\theta }_j+\mathrm{\Delta }{\theta }_j)-y_i(t,{\theta }_j)}{y_i(t,{\theta }_j)}/\phantom{\frac{y_i(t,{\theta }_j+\mathrm{\Delta }{\theta }_j)-y_i(t,{\theta }_j)}{y_i(t,{\theta }_j)}\frac{\mathrm{\Delta }{\theta }_j}{{\theta }_j}}\frac{\mathrm{\Delta }{\theta }_j}{{\theta }_j}.

The time-averaged sensitivity coefficient can be calculated as follows,

(18)S_{ij}=\frac{1}{T} {\int }_0^T \left|S_{ij}(t)\right| \mathrm{d}t\approx \frac{1}{K} \sum _{k=1}^K\left|S_{ij}(t_k)\right|,

where \left\{t_k,k=1,\dots ,K\right\} is an equidistant sub-division of [0,1].

Based on the results of the above parameter sensitivity analysis, the values of key parameters can be changed in a wider range to examine the robustness of the model. For example, for each parameter, we can change its value to 0.1 times, 2.0 times, 5.0 times, 10.0 times, 50.0 times and 100.0 times, respectively, and at the same time keep the other parameters invariant. Then we can check how the system responds to changes in the above parameters. If the model still maintains a small fluctuation for a large range of value changes of a parameter, then it is considered that the model is robust with respect to this parameter.

6 Application of multiscale modeling in bone regeneration and cancer biology

The bone regeneration in tissue engineering is a complex biomedicine problem, which involves scaffolds, cells, signaling molecules, vascularization and their interplay [44]. Many in vivo and in vitro experiments have investigated strategies for promoting bone formation and angiogenesis [39, 44]. However, due to the complexity of bone regeneration within scaffolds, these interwoving biological processes are often investigated in isolation while rarely studied interactively and systematically.

In recent years, many mathematical models of bone regeneration have been established [57]. Geris et al. [16-18] proposed a continuum model that uses a set of PDEs to describe the temporal and spatial evolutions of cell densities and growth factor concentrations, but this model did not include the release of exogenous growth factors. Sanz-Herrera et al. [52-54] employed finite element method to establish a micro-macro numerical model of bone regeneration, integrating the tissue level (i.e., macro level) and the scaffold pore level (i.e., micro level). Checa et al. [8, 9, 51] developed a mechanical-biological model of cell activity by using a lattice method and investigated the dependence of bone regeneration and vascularization on cell implantation and mechanical loading. However, the above models did not consider the process of release of exogenous growth factors from the scaffold nor the response of cells to these growth factors. To overcome these shortcomings, we proposed a continuous-discrete hybrid model (multi-scale PDE-Agent model) that integrated the growth factor release process and the cell response into the vascularized bone regeneration system [61]. This study simulated the response of osteogenesis to growth factor release based on molecular mechanisms and incorporated angiogenesis and nutrient transport processes. It revealed that the combination of BMP2, Wnt and VEGF could better promote angiogenesis and bone formation compared to individual growth factors at the same dose.

Bone remodeling is a long-term physiological process in which bone formation and bone resorption balance each other in the late stage of bone regeneration. It involves stem cell-driven cell lineage development, intercellular interactions between various cell phenotypes and intracellular signaling pathways. In recent years, many mathematical models have been developed to describe bone remodeling. Komarova et al. [32] developed a mathematical model to simulate dynamic changes in cell numbers and bone mass during bone remodeling, which considered the autocrine and paracrine signaling interactions between osteoblasts and osteoclasts in a simplified manner. Lemaire et al. [36] proposed another cell population dynamics model to explain the interaction between osteoblasts and osteoclasts by incorporating the RANK-RANKL-OPG intercellular signaling pathway. Later, Pivonka et al. [48] expanded this signaling pathway and specifically studied the influence of RANK/OPG expression morphology on bone function. In recent years, Pivonka et al. [48, 49] investigated possible treatment options for bone loss caused by imbalance of RANK-RANKL-OPG regulation [49]. However, these models lacked intracellular signaling pathways so they couldnot represent intracellular molecular mechanisms. In addition, these models did not examine the therapeutic effects of cytokines (especially cytokine combinations) on bone remodeling or bone healing. Our study [62] has established a systems model based on intracellular signaling pathways to investigate the roles of cytokines (including Wnt, BMP2 and TGFβ and their combinations) released from scaffolds in controlling the balance of bone formation and resorption. Firstly, a set of ODEs was established to describe the intracellular signaling pathways regulating bone differentiation; then the intracellular signaling pathways and intercellular signaling pathways affecting osteoclast differentiation were coupled, which were integrated into a set of SDEs describing the dynamic changes in the densities of cells belonging to different populations. The synergistic effect of cytokine combinations was quantitatively evaluated using Loewe and Bliss indexes. The results demonstrated that the combination of some growth factors, such as BMP2 and Wnt, and the optimization of their dosage can benefit the effect of bone remodeling.

The above multi-scale modeling approaches can also be applied to cancer biology. Cancer is a very complex and deadly disease. Although biologists have obtained many experimental data at the molecular, cellular, microenvironmental, and tissue levels, few researchers have integrated these data into multiscale models to study tumor growth and its response to drug therapy.

Cellular automata methods have been widely used to simulate tumor growth [24, 47]. Although cellular automata methods have advantages in simulating cell-cell interactions and cell-microenvironment interactions, it’s hard to use this discrete type of approach to describe the reaction-diffusion process of molecules in the tumor microenvironment. While the continuous models that use a set of PDEs to simulate tumor growth [40, 71], however, are difficult to describe interactions between individual cells, for instance the competition for nutrients between cells. Generally speaking, due to the complexity of cancer, neither discrete nor continuous methods alone can accurately simulate the spatio-temporal evolution of cancer. Anderson et al. [3, 50] proposed a hybrid discrete-continuum (HDC) model, which coupled a cellular automaton module with a continuous PDE module. However, their model did not take into account molecular signaling networks. ABM includes gene-protein signaling pathways (such as EGFR signaling pathway) into each cell [68, 69], enabling each cell to select their phenotypic traits based on the state of the signaling pathway [65, 68, 70]. However, previous ABMs did not simulate tumor-induced angiogenesis. We proposed a multi-scale PDE-Agent model [14] to simulate vascularized tumor growth and to study the response of the tumor cells to tyrosine kinase inhibitors by using a set of PDEs to describe the diffusion of VEGF and nutrients in the tumor microenvironment, thus bridging the two interacting processes of tumor growth and angiogenesis. The model was then used to examine angiogenesis-induced drug resistance in brain tumors and predicted a synergistic combination between inhibitors of EGFR and VEGFR [38].

In past decades, some mathematical models of apoptotic regulation have been established. For example, a Boolean model of apoptosis regulation was proposed [55] to quantitatively analyze the internal and external signaling pathways of apoptosis regulation. Based on the law of mass action or Michaelis-Menten kinetics, other continuous models have also been developed for apoptosis pathways. A kinetic model of the signaling pathway controlling apoptosis [35] revealed that inhibition of Caspase 3 and Caspase 9 could lead to a potential positive feedback, resulting in bistability. Later, a mathematical model of the Src-controlled mitochondrial pathway for apoptosis was proposed [5], which was fitted to experimental data and then used to theoretically design optimal treatment strategies. However, these models do not examine the interplays between psychological stress, apoptosis, and drug resistance. Actually, it has long been known that psychological stress and cancers are inextricably linked. For example, it has been reported that psychological stress might promote cancer progression. However, the causal relationship between psychological stress and cancer has been poorly understood, largely due to that the mechanisms by which stress affects cancer development and drug resistance remain largely unclear. To this end, we developed an ODEs network model [58] to describe the dynamics of anti-apoptotic signaling pathways and change in number of cancer cells, which theoretically justified the recent findings that transcription factor Mcl-1 plays an important role in anti-apoptotic regulation. In addition, the model further revealed that the psychological stress of cancer patients affected the switch of synergism patterns of combinations of anticancer drugs.

7 Perspective

We have introduced the two methods of multi-scale modeling (i.e., multiscale SDE model, PDE-Agent model) and demonstrated their applications in bone regeneration and cancer biology. They have good extensibility and can be flexibly applied for modeling of many multi-scale biological systems. On the other hand, the above multi-scale mathematical models still have some limitations and further works are required. For example, the signaling pathways above are mostly based on the knowledge that has been verified by experiments. However, how to infer signaling pathways or networks based on high-throughput data (e.g., proteomics data or transcriptomics data) is an important and basic question [10, 60, 63, 66, 67]. Furthermore, the biological scales in the above models mainly include molecular scale, cellular scale, microenvironmental scale and tissue level, while the population level of clinical and epidemiological concern has not been involved. How to bridge the gap between biological mechanisms and patient population is an important and challenging question for mathematical modeling. We are developing multi-scale stochastic models connecting molecular/cellular mechanisms to population dynamics, using SDEs or PDEs and Monte Carlo simulations to predict the probability distribution of survival time of cancer patients in an in-silico cohort [59, 72].

Acknowledgements

This work was supported by the National Natural Science Foundation of China (Grant Nos. 62273364, 11871070 & 11371060) and the Guangdong Basic and Applied Basic Research Foundation (2020B1515020047).

References

1 Aldridge B B, Burke J M, Lauffenburger D A. . Physicochemical modelling of cell signalling pathways. Nat Cell Biol 2006; 8(11): 1195–1203

2 AlonU. An Introduction to Systems Biology: Design Principles of Biological Circuits. Boca Raton: Chapman & Hall/CRC, 2006

3 Anderson A R A. A hybrid mathematical model of solid tumour invasion: the importance of cell adhesion. Math Med Biol J IMA 2005; 22(2): 163–186

4 Anderson A R A, Chaplain M A J. Continuous and discrete mathematical models of tumor-induced angiogenesis. Bull Math Biol 1998; 60(5): 857–899

5 Ballesta A, Lopez J, Popgeorgiev N. . Data-driven modeling of SRC control on the mitochondrial pathway of apoptosis: implication for anticancer therapy optimization. PLoS Comput Biol 2013; 9(4): e1003011

6 Bonabeau E. Agent-based modeling: Methods and techniques for simulating human systems. Proc Nat Acad Sci USA 2002; 99(Suppl 3): 7280–7287

7 ChapmanS J. Multiscale mathematical modelling in medicine and biology. In: Proceedings of the 18th World IMACS Congress and MODSIM09 International Congress on Modelling and Simulation: Interfacing Modeling and Simulation with Mathematcial and Computational Sciences, Anderssen R S, Bradock R D and Newham L T H eds. Cairns, Australia, 2009, 13–22

8 Checa S, Prendergast P J. A mechanobiological model for tissue differentiation that includes angiogenesis: a lattice-based modeling approach. Ann Biomed Eng 2009; 37(1): 129–145

9 Checa S, Prendergast P J. Effect of cell seeding and mechanical loading on vascularization and tissue formation inside a scaffold: a mechano-biological model using a lattice approach to simulate cell activity. J Biomech 2010; 43(5): 961–968

10 Cheng J, Zhang J, Wu Z. . Inferring microenvironmental regulation of gene expression from single-cell RNA sequencing data using scMLnet with an application to COVID-19. Briefings Bioinf 2021; 22(2): 988–1005

11 ChopardBDroz M. Cellular Automata Modeling of Physical Systems. Cambridge: Cambridge University Press, 1998

12 Chou T C. Derivation and properties of Michaelis-Menten type and Hill type equations for reference ligands. J Theor Biol 1976; 59(2): 253–276

13 Dada J O, Mendes P. Multi-scale modelling and simulation in systems biology. Inter Biol 2011; 3(2): 86–96

14 Ellis P M, Morzycki W, Melosky B. . The role of the epidermal growth factor receptor tyrosine kinase inhibitors as therapy for advanced, metastatic, and recurrent non-small-cell lung cancer: a Canadian national consensus statement. Curr Oncol 2009; 16(1): 27–48

15 Garfinkel A, Tintut Y, Petrasek D. . Pattern formation by vascular mesenchymal cells. Proc Nat Acad Sci USA 2004; 101(25): 9247–9250

16 Geris L, Reed A A C, Vander Sloten J. . Occurrence and treatment of bone atrophic non-unions investigated by an integrative approach. PLoS Comput Biol 2010; 6(9): e1000915

17 Geris L, Vander Sloten J, Van Oosterwyck H. In silico biology of bone modelling and remodelling: regeneration. Phil Trans Royal Soc A: Math Phys Eng Sci 2009; 367(1895): 2031–2053

18 Geris L, Vander Sloten J, Van Oosterwyck H. Connecting biology and mechanics in fracture healing: an integrated mathematical modeling framework for the study of nonunions. Biomech Model Mechanobiol 2010; 9(6): 713–724

19 Gierer A, Meinhardt H. A theory of biological pattern formation. Kybernetik 1972; 12(1): 30–39

20 Gillespie D T. Exact stochastic simulation of coupled chemical reactions. J Phys Chem 1977; 81(25): 2340–2361

21 Gillespie D T. Stochastic simulation of chemical kinetics. Annu Rev Phys Chem 2007; 58(1): 35–55

22 Gordeladze J O, Reseland J E, Duroux-Richard I. . From stem cells to bone: phenotype acquisition, stabilization, and tissue engineering in animal models. ILAR J 2009; 51(1): 42–61

23 Higham D J. An algorithmic introduction to numerical simulation of stochastic differential equations. SIAM Rev 2001; 43(3): 525–546

24 Kansal A R, Torquato S, Harsh IV G R. . Cellular automaton of idealized brain tumor growth dynamics. Biosystems 2000; 55(1/2/3): 119–127

25 Karplus M, Petsko G A. Molecular dynamics simulations in biology. Nature 1990; 347(6294): 631–639

26 Katare S, Bhan A, Caruthers J M. . A hybrid genetic algorithm for efficient parameter estimation of large kinetic models. Comp Chem Eng 2004; 28(12): 2569–2581

27 Kauffman S A. Metabolic stability and epigenesis in randomly constructed genetic nets. J Theor Biol 1969; 22(3): 437–467

28 KitanoH. Foundations of Systems Biology. Cambridge, MA: MIT Press, 2001

29 KlebanerF C. Introduction to Stochastic Calculus with Applications, 2nd ed. Singapore: Imperial College Press, 2005

30 Klepeis J L, Lindorff-Larsen K, Dror R O. . Long-timescale molecular dynamics simulations of protein structure and function. Curr Opin Struct Biol 2009; 19(2): 120–127

31 KlippELiebermeister WWierlingC, . Systems Biology: A Textbook, 1st ed. Weinheim: Wiley-VCH, 2009, 13–22

32 Komarova S V, Smith R J, Dixon S J. . Mathematical model predicts a critical role for osteoclast autocrine regulation in the control of bone remodeling. Bone 2003; 33(2): 206–215

33 Komori T. Regulation of bone development and maintenance by Runx2. Front Biosci 2008; 13(3): 898–903

34 Lee J M, Gianchandani E P, Eddy J A. . Dynamic analysis of integrated signaling, metabolic, and regulatory networks. PLoS Comput Biol 2008; 4(5): e1000086

35 Legewie S, Blüthgen N, Herzel H. Mathematical modeling identifies inhibitors of apoptosis as mediators of positive feedback and bistability. PLoS Comput Biol 2006; 2(9): e120

36 Lemaire V, Tobin F L, Greller L D. . Modeling the interactions between osteoblast and osteoclast activities in bone remodeling. J Theor Biol 2004; 229(3): 293–309

37 Liang J, Qian H. Computational cellular dynamics based on the chemical master equation: a challenge for understanding complexity. J Comput Sci Tech 2010; 25(1): 154–168

38 Liang W, Zheng Y, Zhang J. . Multiscale modeling reveals angiogenesis-induced drug resistance in brain tumors and predicts a synergistic drug combination targeting EGFR and VEGFR pathways. BMC Bioinf 2019; 20(Suppl 7): 203

39 Lovett M, Lee K, Edwards A. . Vascularization strategies for tissue engineering. Tissue Eng Part B: Reviews 2009; 15(3): 353–370

40 Macklin P, McDougall S, Anderson A R. . Multiscale modelling and nonlinear simulation of vascular tumour growth. J Math Biol 2009; 58(4/5): 765–798

41 Michaelis L, Menten M L. Die kinetik der invertinwirkung. Biochem Z 1913; 49: 333–369

42 MortonK WMayers D F. Numerical Solution of Partial Differential Equations: An Introduction. Cambridge: Cambridge University Press, 2005

43 Neumann L, Pforr C, Beaudouin J. . Dynamics within the CD95 death-inducing signaling complex decide life and death of cells. Mol Syst Biol 2010; 6(1):

44 Nguyen L H, Annabi N, Nikkah M. . Vascularized bone tissue engineering: approaches for potential improvement. Tissue Eng Part B: Reviews 2009; 18(5): 363–382

45 NicolisGPrigogine I. Self-organization in Nonequilibrium Systems: From Dissipative Structures to Order Through Fluctuations. New York: John Wiley & Sons, 1977

46 ØksendalB K. Stochastic Differential Equations: An Introduction with Applications, 5th ed. New York: Springer-Verlag, 2003, 61–72

47 Patel A A, Gawlinski E T, Lemieux S K. . A cellular automaton model of early tumor growth and invasion: the effects of native tissue vascularity and increased anaerobic tumor metabolism. J Theor Biol 2001; 213(3): 315–331

48 Pivonka P, Zimak J, Smith D W. . Model structure and control of bone remodeling: a theoretical study. Bone 2008; 43(2): 249–263

49 Pivonka P, Zimak J, Smith D W. . Theoretical investigation of the role of the RANK-RANKL-OPG system in bone remodeling. J Theor Biol 2010; 262(2): 306–316

50 Rejniak K A, Anderson A R A. Hybrid models of tumor growth. Wiley Interdiscip Rev: Syst Biol Med 2011; 3(1): 115–125

51 Sandino C, Checa S, Prendergast P J. . Simulation of angiogenesis and cell differentiation in a CaP scaffold subjected to compressive strains using a lattice modeling approach. Biomaterials 2009; 31(8): 2446–2452

52 Sanz-Herrera J A, Garcia-Aznar J, Doblare M. Micro-macro numerical modelling of bone regeneration in tissue engineering. Comput Meth Appl Mech Eng 2008; 197(33): 3092–3107

53 Sanz-Herrera J A, Garcia-Aznar J, Doblare M. Simulation of bone remodelling and bone ingrowth within scaffolds. Key Eng Mater 2008; 377: 225–273

54 Sanz-Herrera J A, García-Aznar J M, Doblaré M. A mathematical approach to bone tissue engineering. Phi Trans Roy Soc A: Math Phys Eng Sci 2009; 367(1895): 2055–2078

55 Schlatter R, Schmich K, Avalos Vizcarra I. . ON/OFF and beyond-a Boolean model of apoptosis. PLoS Comput Biol 2009; 5(12): e1000595

56 Stokes C L, Lauffenburger D A. Analysis of the roles of microvessel endothelial cell random motility and chemotaxis in angiogenesis. J Theor Biol 1991; 152(3): 377–403

57 Sun W, Lal P. Recent development on computer aided tissue engineering—a review. Comput Meth Prog Biomed 2002; 67(2): 85–103

58 Sun X Q, Bao J G, Nelson K C. . Systems modeling of anti-apoptotic pathways in prostate cancer: psychological stress triggers a synergism pattern switch in drug combination therapy. PLoS Comput Biol 2013; 9(12): e1003358

59 Sun X Q, Bao J G, Shao Y Z. Mathematical modeling of therapy-induced cancer drug resistance: connecting cancer mechanisms to population survival rates. Sci Rep 2016; 6(1): 22498

60 Sun X Q, Hu B. Mathematical modeling and computational prediction of cancer drug resistance. Briefings Bioinf 2018; 19(6): 1382–1399

61 Sun X Q, Kang Y Q, Bao J G. . Modeling vascularized bone regeneration within a porous biodegradable CaP scaffold loaded with growth factors. Biomaterials 2013; 34(21): 4971–4981

62 Sun X Q, Su J, Bao J G. . Cytokine combination therapy prediction for bone remodeling in tissue engineering based on the intracellular signaling pathway. Biomaterials 2012; 33(33): 8265–8276

63 Sun X Q, Zhang J, Nie Q. Inferring latent temporal progression and regulatory networks from cross-sectional transcriptomic data of cancer samples. PLoS Comput Biol 2021; 17(3): e1008379

64 Sun X Q, Zhang J, Zhao Q. . Stochastic modeling suggests that noise reduces differentiation efficiency by inducing a heterogeneous drug response in glioma differentiation therapy. BMC Syst Biol 2016; 10(1): 73

65 Wang Z H, Zhang L, Sagotsky J. . Simulating non-small cell lung cancer with a multiscale agent-based model. Theor Biol Med Model 2007; 4: 50

66 Zhang J, Guan M G, Wang Q L. . Single-cell transcriptome-based multilayer network biomarker for predicting prognosis and therapeutic response of gliomas. Briefings Bioinf 2020; 21(3): 1080–1097

67 Zhang J J, Zhu W B, Wang Q L. . Differential regulatory network-based quantification and prioritization of key genes underlying cancer drug resistance based on time-course RNA-seq data. PLoS Comput Biol 2019; 15(11): e1007435

68 Zhang L, Athale C A, Deisboeck T S. Development of a three-dimensional multiscale agent-based tumor model: simulating gene-protein interaction profiles, cell phenotypes and multicellular patterns in brain cancer. J Theor Biol 2007; 244(1): 96–107

69 Zhang L, Chen L L, Deisboeck T S. Multi-scale, multi-resolution brain cancer modeling. Mathe Comp Simul 2009; 79(7): 2021–2035

70 Zhang L, Wang Z H, Sagotsky J A. . Multiscale agent-based cancer modeling. J Math Biol 2009; 58(4/5): 545–559

71 Zheng X, Wise S M, Cristini V. Nonlinear simulation of tumor necrosis, neo-vascularization and tissue invasion via an adaptive finite-element/level-set method. Bull Math Biol 2005; 67(2): 211–259

72 Zheng Y J, Bao J G, Zhao Q Y. . A spatio-temporal model of macrophage-mediated drug resistance in glioma immunotherapy. Mol Cancer Ther 2018; 17(4): 814–824