By Luis L. Bonilla
The collective motion of cells is very important for physiological processes—such as organ growth and regeneration, tissue repair, and wound healing—as well as pathological processes like cancer and retinopathies. Here, we will first explain mathematical models of angiogenesis, or the multiscale process of blood vessel growth from arteries and veins. We will then address the motion of epithelial cells on a planar substrate, which contributes to our understanding of morphogenesis and the invasion of normal tissue by cancerous cells.
Angiogenesis is triggered by hypoxia, a lack of oxygen for cells in some area of tissue. Hypoxic cells secrete substances called vessel endothelial growth factors (VEGFs) that diffuse and reach a primary blood vessel. The vessel opens its wall in response, and endothelial cells emerge to form a capillary sprout of cells that collectively move towards the VEGF source. The sprout eventually becomes a new vessel that carries blood with oxygen and nutrients toward the hypoxic cells. Once the blood and oxygen reach the hypoxic region, it stops secreting growth factors and may begin to secrete anti-angiogenic substances. A regular vessel network may then form after capillaries with insufficient blood flow are pruned. Under normal physiological conditions, angiogenic and anti-angiogenic activities balance each other out [6, 10]; imbalance often results in disease, including cancer. In fact, after a tumor within a tissue reaches two millimeters in size, it requires additional nutrients and oxygen to continue its growth. Its hypoxic cells thus secrete growth factors and induce angiogenesis. Unlike normal cells, cancerous cells continue to issue growth factors and attract blood vessels, which also supply them with a handy transportation system to other organs in the body.
Figure 1. Hybrid stochastic model of tumor-induced angiogenesis. 1a. A network of blood vessels simulated by the stochastic model. The level curves of the density of the vessel endothelial growth factor are depicted by the red lines. 1b. Comparison of the active tip density profile \(\tilde{p}(t,x,0)\) (calculated by averaging the stochastic description over 400 replicas) to the profile of the moving soliton. Figure adapted from [5].
The complex angiogenesis process involves serial changes across multiple scales, such as communication between cells via their membranes; differentiation of the roles of cells at the tip of the sprout versus those in the wake of the tip (respectively called tip and stalk cells); sprout branching, fusion, and pruning; formation of sprout walls; and circulation of blood through the new vessels.
The study of angiogenesis requires a combination of experiments, theories, and mathematical and computational modeling. Mathematical models that target specific aspects at subcellular, cellular, or mesoscopic scales [4, 7, 8] may involve continuum tools (e.g., partial differential equations (PDEs) for cellular densities, elastic and phase fields, and the diffusion and reaction of substances and growth factors) and discrete tools (e.g., agent-based models or cellular Potts models (CPMs) for subcellular and cellular mechanisms or extracellular space). Hybrid models combine continuum and discrete approaches to encompass diverse length and time scales; they often involve stochastic differential equations and stochastic birth and death processes.
Animation 1. Progression of wet age-related macular degeneration due to capillaries that issue from choroid blood vessels and are attracted to hypoxic cells in Bruch’s membrane (which is deformed by lipid deposits) and hypoxic photoreceptors at \(y=L\). Animation courtesy of [12].
While discrete and hybrid models enable the incorporation of detailed biological mechanisms, continuum approaches are more amenable to analysis and numerical simulation. Averaging discrete or hybrid models over small scales may produce continuum equations, which one can further analyze and relate back to the original model. As an example, consider the numerical simulation of the two-dimensional (2D) hybrid model in Figure 1a. Stalk cells simply follow the tip cells (mass particles) to build each growing blood capillary vessel, which coincides with the tip cells’ trajectory. At random times, new tip cells branch out of existing ones with some probability. When a particle collides with another’s trajectory, it stops and disappears; this corresponds to anastomosis (the fusion of vessels). The position \(\mathbf{X}^i\) and velocity \(\mathbf{v}^i\) of particle \(i\) satisfy the Ito differential equations:
\[d\mathbf{X}^i(t) = \mathbf{v}^i(t)\,dt, \\ d\mathbf{v}^i(t) = \beta\left[-\mathbf{v}^i(t) + \mathbf{F}\!\left(C(t,\mathbf{X}^i(t))\right)\! \,\,\right] dt+ \sqrt{\beta}\, d\mathbf{W}^i(t). \tag1\]
The forces that act on particle \(i\) are friction \(-\beta\mathbf{v}^i\), a chemotactic force \(\mathbf{F}\propto\nabla C\), and a random force that is proportional to the differential of a Wiener process \(\mathbf{W}^i(t)\). The number of particles at time \(t\) is \(N(t)\), and the VEGF density \(C(t,\mathbf{x})\) diffuses and is consumed by the moving particles. A number of active tips initially emerge from the \(y\) axis and move towards the tumor that is releasing VEGFs. Numerical simulations show that in the simple slab geometry of Figure 1, the density of active particles advances as a solitary wave towards the tumor. Averaging over a sufficiently large number of replicas of the stochastic process, the density \(p(t,\mathbf{x},\mathbf{v})\) of active cells at \(\mathbf{x}\) that move with velocity \(\mathbf{v}\) obeys the following PDE:
\[\frac{\partial p}{\partial t}(t,\mathbf{x},\mathbf{v})= \alpha(C(t,\mathbf{x}))\,
p(t,\mathbf{x},\mathbf{v})\delta_{\sigma_v}(\mathbf{v}-\mathbf{v}_0) - \Gamma\, p(t,\mathbf{x},\mathbf{v}) \int_0^t \tilde{p}(s,\mathbf{x})\, ds - \\
\mathbf{v}\cdot\nabla_x p(t,\mathbf{x},\mathbf{v}) - \beta\nabla_v \cdot \{[\mathbf{F}(C(t,\mathbf{x}))-\mathbf{v}]\, p(t,\mathbf{x},\mathbf{v})\} + \frac{\beta}{2} \Delta_{v} p(t,\mathbf{x},\mathbf{v}). \tag2\]
The integral of \(p\) over velocity is the number density of active tip cells \(\tilde{p}(t,\mathbf{x})\). The first two terms on the righthand side of \((2)\) respectively describe vessel tip branching and anastomosis; the following components are the convective and diffusive terms that are associated with the Fokker-Planck equation for the probability density that corresponds to \((1)\). Under the slab geometry of Figure 1, the density of active tips \(\tilde{p}(t,\mathbf{x})\) volves towards a solitary wave that moves on the horizontal axis \(y=0\). A dominant balance of time rate, convection, and source terms in \((2)\) produces
\[\tilde{p}_s(t,x)=\frac{(2K\Gamma+\mu^2)c}{2\Gamma(c-F_x)}\mbox{sech}^2\!\left[\frac{\sqrt{2K\Gamma+\mu^2}}{2(c-F_x)}(x-X(t))\right]\!,\quad \dot{X}\equiv\frac{dX}{dt}=c. \tag3\]
This result is the soliton wave of the Korteweg-de Vries equation for water waves with velocity \(c\) and an additional shape parameter \(K\). These parameters are called collective coordinates and slowly evolve according to ordinary differential equations (ODEs) that reflect the variation of VEGF concentration and the subdominant terms in \((2)\). \(F_x\) is the \(x\) component of the chemotactic force and \(\mu\) is a function of the branching function \(\alpha\). Figure 1b compares the average over 400 realizations of the stochastic process to the moving soliton in \((3)\) with slowly varying shape and velocity. The soliton effectively approximates both the solutions of the PDE in \((2)\) and the ensemble average over replicas of the stochastic process after a formation stage and before it reaches the VEGF source. The soliton is thus an attractor for this intermediate stage. The center of mass \(X(t)\) is a good measure of the advancing angiogenic network, and its collective coordinates—which provide a very simple description in terms of ODEs—are a suitable starting point for the analysis of the control of angiogenesis in simple models.
Figure 2. Hybrid discrete models of cellular motion. 2a. Three successive snapshots of the progression of wet age-related macular degeneration due to capillaries that issue from choroid blood vessels and are attracted to hypoxic cells in Bruch’s membrane (which is deformed by lipid deposits) and hypoxic photoreceptors at \(y=L\). Pink represents tip cells and red represents stalk cells. Vessel endothelial growth factor level curves are depicted with a color gradient. See Animation 1 for a video of the numerical simulation. 2b. A population of cancerous epithelial cells (in magenta) invades healthy tissue (in green) and creates a rough interface with pockets of other cells. See Animation 2 for a video of the numerical simulation. Figure 2a adapted from [12] and 2b adapted from [3].
Cellular models of angiogenesis combine descriptions of cells and their interactions with the continuum fields that are responsible for their motion, including gradients of chemicals (chemotaxis), adhesion (haptotaxis), and elastic strain (durotaxis). CPM is a widely used and flexible but computationally expensive approach to this problem. It consists of a square or cubic grid whose elementary units \(\sigma\) (sites) have a label \(\tau\) that characterizes them as endothelial cells, extracellular matrices, fluids, or so forth. An energy functional determines the overall configuration; a 2D example is
\[H = \sum_{\rm sites} J_{\tau,\tau'} (1-\delta_{\sigma\sigma'})+ \sum_{\rm cells}\gamma_\tau (a_\sigma-A_\sigma)^2 + \sum_{\rm cells} \sum_{\rm sites} \rho_\tau (p_\sigma-P_\sigma)^2+ H_\text{chem}.\tag4\]
The configuration follows Monte Carlo dynamics. At discrete times, a random site is selected to change; the difference between the new and previous configurations \(\Delta H\) is accepted if negative and accepted with probability \(e^{-\Delta H}\) if positive. All moduli \(J_{\tau,\tau'}\), \(\gamma_\tau\), and \(\rho_\tau\) are positive. The first term in \((4)\) facilitates adhesion between sites of the same cell type. The area \(a_\sigma\) and perimeter \(p_\sigma\) of a cell attempt to reach target values \(A_\sigma\) and \(P_\sigma\), respectively. During chemotaxis, \(\Delta H_{\rm chem}=-\mu_\sigma[C(t,\mathbf{x})-C(t,\mathbf{x}')]\); here, \(\mathbf{x}\) and \(\mathbf{x}'\) are two randomly selected neighboring lattice cells and \(\mu_\sigma>0\) is the chemical potential. As in the simpler mesoscopic model, the VEGF concentration \(C(t,\mathbf{x})\) obeys a reaction-diffusion equation. Including additional constraints allows us to restrict the length of the cells, incorporate the equations of elasticity, and add a durotaxis term in \((4)\) to consider the tensions that act on cells as they move on a substrate. To distinguish between tip and stalk cells, we can add ODEs for each cell that mimic Notch signaling: the way in which the cells react to the external VEGFs based on proteins and receptor activation. Some cells with a high concentration of a certain protein become tip cells, exhibit more sensitivity to external fields, and can move out of an existing capillary if the external strain fields allow [13].
Animation 2. A population of cancerous epithelial cells (in magenta) invades healthy tissue (in green) and creates a rough interface with pockets of other cells. Animation courtesy of [3].
Figure 2a uses these ideas to simulate the wet phase of age-related macular degeneration (AMD): the leading cause of blindness among people over 60 years old. In general, the photoreceptors in the macula that are responsible for central vision are separated from the choroid blood vessels by a membrane and retinal pigment epithelium cells. This membrane gradually becomes more rigid and thick with age, losing its ability to filter oxygen and nutrients from the choroid and instead allowing debris to deposit on it. The hypoxic photoreceptors in patients with wet AMD release extra VEGFs that induce angiogenesis from the choroid. The new blood vessels may cross the membrane, leak blood, and knock down photoreceptors — thereby provoking the loss of central vision.
Vertex models are less computationally expensive than CPMs and are widely applied to study the motion of epithelial cells. In these models, the polygonal cells that tile the plane with a Voronoi tessellation have an energy functional without adhesion or chemotaxis. The dual Delaunay triangulation consists of cell centers whose motion can manifest in a variety of ways within the self-propelled Voronoi vertex models [1, 2]. Based on the selection of parameters in the energy functional, cellular tissue can behave as an elastic solid or a fluid. Figure 2b depicts the invasion of a monolayer tissue of healthy cells by similar but cancerous cells [3].
Several successful continuum models of collective cell migration in tissues describe supracellular length and time scales and are more amenable to analysis than their subcellular and cellular counterparts. Their dynamics are governed by PDEs for fields such as average density, velocity, polarity, and nematic order parameter [1]. Deriving these continuum models from subcellular and cellular models—e.g., CPMs or active vertex models—is an important but unexplored task. Coarse graining procedures that utilize Poisson brackets seem promising [9, 11], but much more work is needed.
Luis Bonilla delivered a minisymposium presentation on this research at the 2023 SIAM Conference on Computational Science and Engineering, which took place in Amsterdam, the Netherlands, earlier this year.
References
[1] Alert, R., & Trepat, X. (2020). Physical models of collective cell migration. Ann. Rev. Cond. Matt. Phys., 11, 77-101.
[2] Barton, D.L., Henkes, S., Weijer, C.J., & Sknepnek, R. (2017). Active vertex model for cell-resolution description of epithelial tissue mechanics. PLoS Comput. Biol., 13(6), e1005569.
[3] Bonilla, L.L., Carpio, A., & Trenado, C. (2020). Tracking collective cell motion by topological data analysis. PLoS Comput. Biol., 16(12), e1008407.
[4] Bonilla, L.L., Carretero, M., & Terragni, F. (2019). Stochastic models of blood vessel growth. In IHPStochDyn 2017: Stochastic dynamics out of equilibrium (pp. 413-436). Paris, France: Springer.
[5] Bonilla, L.L., Carretero, M., Terragni, F., & Birnir, B. (2016). Soliton driven angiogenesis. Sci. Rep., 6, 31296.
[6] Carmeliet, P. (2005). Angiogenesis in life, disease and medicine. Nature, 438(7070), 932-936.
[7] Flegg, J.A., Menon, S.N., Byrne, H.M., & McElwain, D.L.S. (2020). A current perspective on wound healing and tumour-induced angiogenesis. Bull. Math. Biol., 82(2), 23.
[8] Heck, T.A.M., Vaeyens, M.M., & Van Oosterwyck, H. (2015). Computational models of sprouting angiogenesis and cell migration: Towards multiscale mechanochemical models of angiogenesis. Math. Model. Nat. Phenom., 10(1), 108-141.
[9] Hernandez, A., & Marchetti, M.C. (2021). Poisson-bracket formulation of dynamics of fluids of deformable particles. Phys. Rev. E, 103(3), 032612.
[10] Szymborska, A., & Gerhardt, H. (2018). Hold me, but not too tight – endothelial cell-cell junctions in angiogenesis. Cold Spring Harb. Perspect. Biol., 10(8), a029223.
[11] Triguero-Platero, G., Ziebert, F., & Bonilla, L.L. (2023). Coarse-graining the vertex model and its response to shear. Preprint, arXiv:2302.04111.
[12] Vega, R., Carretero, M., & Bonilla, L.L. (2021). Anomalous angiogenesis in retina. Biomedicines, 9(2), 224.
[13] Vega, R., Carretero, M., Travasso, R.D.M., & Bonilla, L.L. (2020). Notch signaling and taxis mechanisms regulate early stage angiogenesis: A mathematical and computational model. PLoS Comput. Biol., 16(1), e1006919.
Luis L. Bonilla is a professor of applied mathematics and director of the Gregorio Millán Barbany Institute of Fluid Dynamics, Nanoscience and Industrial Mathematics at Universidad Carlos III de Madrid in Spain.