Multicellular model of angiogenesis

Takashi Nakazawa; Sohei Tasaki; Kiyohiko Nakai; Takashi Suzuki; Takashi Nakazawa; Sohei Tasaki; Kiyohiko Nakai; Takashi Suzuki

doi:10.3934/bioeng.2022004

AIMS Bioengineering

2022, Volume 9, Issue 1: 44-60. doi: 10.3934/bioeng.2022004

Previous Article Next Article

Research article Special Issues

Multicellular model of angiogenesis

1.
Center for Mathematical Modelling and Data Science, Osaka University, 1-3, Machikaneyama-cho, Toyonaka-shi, Osaka 560-8531, Japan
2.
Department of Mathematics, Faculty of Science, Hokkaido University, Kita-10, Nishi-8, Kita-ku, Sapporo 060-0810, Japan
3.
Chugai Pharmaceutical Co. Ltd., 2-1-1 Nihonbashi-Muromachi, Chuo-ku, Tokyo 103-8324, Japan

Received: 12 December 2022 Revised: 23 January 2022 Accepted: 25 January 2022 Published: 28 January 2022

This paper presents a mathematical model governing the dynamics of a morphogenetic vascular endothelial cell (EC) during angiogenesis, and vascular growth formed by EC. Especially, we adopt a multiparticle system for modeling these cells. This model does not distinguish a tip cell from a stalk cell. A formed vessel is modeled using phase-field equation to prevent capillary expansion with time stepping in particular. Numerical simulation reveals that all cells are moving in the direction of high concentration of vascular endothelial growth factor (VEGF), and that they are mutually repellent in cases in which they are closer than some threshold.

Keywords:

Citation: Takashi Nakazawa, Sohei Tasaki, Kiyohiko Nakai, Takashi Suzuki. Multicellular model of angiogenesis[J]. AIMS Bioengineering, 2022, 9(1): 44-60. doi: 10.3934/bioeng.2022004

Related Papers:

[1]	Annabelle Collin, Hadrien Bruhier, Jelena Kolosnjaj, Muriel Golzio, Marie-Pierre Rols, Clair Poignard . Spatial mechanistic modeling for prediction of 3D multicellular spheroids behavior upon exposure to high intensity pulsed electric fields. AIMS Bioengineering, 2022, 9(2): 102-122. doi: 10.3934/bioeng.2022009
[2]	Noorehan Yaacob, Sharidan Shafie, Takashi Suzuki, Mohd Ariff Admon . Signal transduction from ligand-receptor binding associated with the formation of invadopodia in an invasive cancer cell. AIMS Bioengineering, 2022, 9(3): 252-265. doi: 10.3934/bioeng.2022017
[3]	Daniel N. Riahi, Saulo Orizaga . Modeling and computation for unsteady blood flow and solute concentration in a constricted porous artery. AIMS Bioengineering, 2023, 10(1): 67-88. doi: 10.3934/bioeng.2023007
[4]	Sunil Gautam, Sangeeta Kanakraj, Azriel Henry . Computational approach using machine learning modelling for optimization of transesterification process for linseed biodiesel production. AIMS Bioengineering, 2022, 9(4): 319-336. doi: 10.3934/bioeng.2022023
[5]	Sébastien Bernacchi, Simon Rittmann, Arne H. Seifert, Alexander Krajete, Christoph Herwig . Experimental methods for screening parameters influencing the growth to product yield (Y_(x/CH4)) of a biological methane production (BMP) process performed with Methanothermobacter marburgensis. AIMS Bioengineering, 2014, 1(2): 72-87. doi: 10.3934/bioeng.2014.2.72
[6]	Sébastien Bernacchi, Michaela Weissgram, Walter Wukovits, Christoph Herwig . Process efficiency simulation for key process parameters in biological methanogenesis. AIMS Bioengineering, 2014, 1(1): 53-71. doi: 10.3934/bioeng.2014.1.53
[7]	Urooj Ainuddin, Maria Waqas . Finite state machine and Markovian equivalents of the lac Operon in E. coli bacterium. AIMS Bioengineering, 2022, 9(4): 400-419. doi: 10.3934/bioeng.2022029
[8]	Maria Waqas, Urooj Ainuddin, Umar Iftikhar . An analog electronic circuit model for cAMP-dependent pathway—towards creation of Silicon life. AIMS Bioengineering, 2022, 9(2): 145-162. doi: 10.3934/bioeng.2022011
[9]	Sarbaz H. A. Khoshnaw, Kawther Y. Abdulrahman, Arkan N. Mustafa . Identifying key critical model parameters in spreading of COVID-19 pandemic. AIMS Bioengineering, 2022, 9(2): 163-177. doi: 10.3934/bioeng.2022012
[10]	Ayub Ahmed, Bashdar Salam, Mahmud Mohammad, Ali Akgül, Sarbaz H. A. Khoshnaw . Analysis coronavirus disease (COVID-19) model using numerical approaches and logistic model. AIMS Bioengineering, 2020, 7(3): 130-146. doi: 10.3934/bioeng.2020013

Abstract

1. Introduction

Angiogenesis is a multicellular phenomenon by which new blood vessels emerge from an existing vascular system. Especially, tip cell selection plays an important role not only for lateral inhibitation but also for elongation, branching, anastomosis, vessel stabilization, lumen formation, and so on [1], [2]. These events have elucidated the roles of endothelial cell (EC) signaling with the vascular endothelial growth factor (VEGF)-receptor, angiopoitin-Tie2 and Eprin–Eph pathways [1], [3], [4]. Moreover, it has been increasingly clarified that interaction between ECs and mural cells, including vascular smooth muscle cells and pericytes, has been implicated in the maintenance of the angiogenic process for appropriate organization [3]–[7].

In fact, mathematical modeling for vascular growth has attracted considerable interest in the last few years along with elucidation of the mechanisms of development and proliferative and quiescent features, not only of normal cells but also of tumor cells. As a result, an enormous number of phenomena in biochemistry can be subject to mathematical modeling, such as proteins, capillaries, and tumors. Prime examples include the work of Frieboes et al., who coupled phase fields to discrete tumor growth with discrete random walks to model angiogenesis [8]–[10]. Also, Lima et al. [11] proposed the ten-species model. Perfahl et al. [12] considered multiscale modeling of vascular tumor growth.

Cellular automata can be a type of mathematical model that is useful for simulating systems and diseases at the cellular level. Features of cellular automata are simple and discrete systems in which a limited number of states are defined for one cell. Cellular automata can generate complicated behavior based on simple and local rules that determine the states of cells and neighboring cells, whereas distant cells exert no influence [13]. To date, using cellular automata, many models have been presented to elucidate cancer behavior, especially tumor growth and its effective factors [12], [14]–[16]. Competition for nutrients among normal and cancerous cells has been analyzed, as has the influence of the immune system response on tumor growth, as described in an earlier study [17].

Despite the numerous earlier studies, how the spatiotemporal regulation of molecules affects morphogenetic cell movement, and what type of cell collective/cell movement is involved in angiogenesis have remained elusive, largely because of the lack of a stable methodology for visualizing and assessing EC movement during angiogenesis. To clarify the relation between individual cell movements and angiogenesis morphology, and to dissect the underlying molecular and cellular mechanisms, Arima et al. [18] first established a system in which dynamic cell behavior is visualized using time-lapse microscopy. They set out to identify patterns of cellular behavior in an angiogenesis model through computational data processing. As a result, cell movements were quite dynamic. The ECs moved while changing their mutual relative positions at the tip in elongating branches. On the other hand, Takubo et al. [19] analyzed EC behaviours in an in–vitro angiogenic sprouting assay using mouse aortic explants in combination with mathematical modelling. From an experimentally validated mathematical model, cohesive movements with anisotropic cell-to-cell interactions charactarized the EC motolity, which may drive branch elongation depending on a constant cell supply.

In light of the findings described above, we suggest a multilevel model to simulate dynamic cell movement affected by VEGF. Particularly, a multi-agent model is applied to describe cell movement. Some particles with a position vector and polarity mutually interact constantly, with effects of chemotaxis deriving from VEGF. Such cell movement leads to blood vessel formation based on a phase–field equation. The Cahn–Hilliard type, a representative fourth-order partial differential equation, is adopted. As described herein, both mathematical models are coupled to simulate angiogenesis. Then an arbitrary numerical scheme is prepared.

Thermodynamically consistent phase field models refer to a class of models satisfying the second law of thermodynamics. These models have been used for modeling numerous non-equilibrium multi-phase thermodynamical processes ranging from material science and fluid science to life sciences. For this study, we specifically examined a phase field model for time-dependent dynamic of binary material systems in which A is represented by a phase variable at φ = 1 and phase B at φ = 0. We chose the phase variable φ ∈ [0,1] with φ identified at the volume fraction of material A and 1 − φ the volume fraction of material B. The interface between the two phases, known as the transition layer between the two phases in the phase field model, is given as 0 < φ < 1. The interface is defined mathematically at $φ = \frac{1}{2}$ .

The Cahn-Hilliard, introduced by J. Cahn and J. Hilliard in [20] to describe processes of phase separation, includes components of a binary fluid separated and forming pure domains for each component. It can be interpreted as the H⁻¹ gradient flow of the Cahn–Hilliard energy functional as

$F (φ) = \int_{Ω} [\frac{ϵ_{1}}{2} | \nabla φ |^{2} + f (φ)] d x,$ (1.1)

where ϵ₁ is a parameter expressing the strength of the conformational entropy. For immiscible binary materials, we chose the bulk energy density as a double well potential for this study

$f (φ) = ϵ_{2} φ^{2} {(1 - φ)}^{2},$ (1.2)

where ϵ₂ represents the strength of the bulk mixing free energy.

Such a Cahn–Hilliard equation is a typical fourth-order partial differential equation. It is very difficult to obtain a numerical solution stably. A structure-preserving scheme based on the discrete variational derivative method (DVDM) proposed by Furihata and Matsuo [21] has successfully calculated the Cahn–Hilliard equation stably. However, the discontinuous Galerkin method has also been used to trace an interface between material A and material B with high efficiency [22], [23]. Nevertheless, these numerical techniques involve quite cumbersome procedures. Often, we face difficulties in implementing these numerical schemes into our mathematical model with fourth-order partial differential equations. For the study described herein, we chose the Morley finite element method discussed in [24], [25]. The Morley element is implemented in Freefem++ [26] using numerical calculations.

This paper is organized as follows. Section 2 introduces a governing equation representing cell movement as a multi-agent system, with vessel development based on a phase-field model, particularly the Cahn–Hilliard equation. Numerical procedures used to calculate such a coupling model between multi-agent system and phase-field equation are presented in Section 3. Throughout the numerical calculations, numerical results are described in Section 4. Finally, Section 5 concludes this paper.

2. Problem formulation

In this section, we present some basic concepts related to the suggested mathematical multicellular model for angiogenesis. For these purposes, one must describe some mathematical notations.

One can let $Ω \subset ℝ^{d} (d = 1,2,3)$ be a bounded domain, where a position vector is expressed as x ∈ Ω. Cells are located at a position vector x_i and more moving at a velocity vector (polarity) q_i at a time 0 < t ≤ T for a positive constant T, where 0 < i ≤ N_Cell represents an index for distinguishing each cell. In fact, VEGF c attracts such cells to be closer. Consequently, a new blood vessel (capillary) is formed after each cell degenerates fibronectin f and passes. At the next section, a position vector x_i(t), a velocity vector (polarity) q_i(t), and a fibronectin concentration f(t,x) are constructed under a given VEGF c(t,x).

2.1. Governing equation

As discussed above, we suggest equations of three kinds for modeling a mathematical multicellular model for angiogenesis, position vector, polarity for cells, and the fibronectin concentration. Hereinafter, it is assumed that VEFG is given. Also, arbitrary parameters cell mobility a₁, repulsive force between cells and fibronectin a₂ for a position vector x_i and M_q, chemotaxis to VEGF concentration a₃, cell alignment effect a₄, reference size of cell polarity q₀ for polarity q_i and square of transition area width of fibronectin phase variable D, mobility of fibronectin regions M, radius of influence of cells a₅, fibronectin enzyme μ for the fibronectin concentration f are prepared. Let Δt represent a time increment, with n and N respectively denoting a time step for a time t = nΔt and a maximum times step for T = NΔt. These equations are expressed as described hereinafter.

2.1.1. Equation of cell motion

From the discussion presented in the section above, a governing equation of a position vector x_i with the index i distinguishing each cell at a time t, is defined as follows

$\frac{d x_{i}}{d t} = ϕ (f) {{\hat{q}}_{i} + \sum_{j \neq i, j = 1}^{N_{c e l l}} g_{1} (r_{i, j}) \frac{x_{j} - x_{i}}{r_{i, j}}} - a_{2} \nabla f in Ω \times (0, T),$ (2.1)

$x_{i} = x_{i}^{0} a t t = 0.$ (2.2)

The first term on the right-hand side represents that a particle can move in a vessel for f = 0, and not outside of a vessel for f = 1, where ϕ(f) is expressed as

$ϕ (f) = a_{1} e x p (- 4 f),$ (2.3)

and where

${\hat{q}}_{i}$ represents the unit vector of q_i, expressed as

${\hat{q}}_{i} = \frac{q_{i}}{| q_{i} |} .$ (2.4)

The second term on the right-hand side depicts a collision: interaction between particles. Also, r_i_,j, g₁(r_i_,j) studied in [19], [27] are a distance and a restitution coefficient between x_i and x_j as

$r_{i, j} = | x_{j} - x_{i} |,$ (2.5)

$g_{1} (r_{i, j}) = {\begin{array}{l} F_{r e p} \frac{r_{i, j} - d_{c o r e}}{d_{c o r e}} & i f r_{i, j} < d_{c o r e} \\ 0 & i f d_{c o r e} \leq r_{i, j} < d_{n e u t r a l} \\ F_{a t t} \frac{r_{i, j} - d_{n e u t r a l}}{d_{a d h} - d_{n e u t r a l}} & i f d_{n e u t r a l} \leq r_{i, j} < d_{a d h} \\ F_{a t t} \frac{d_{r e a c h} - r_{i, j}}{d_{r e a c h} - d_{a d h}} & i f d_{a d h} \leq r_{i, j} < d_{r e a c h} \\ 0 & i f d_{r e a c h} \leq r_{i, j} \end{array} .$ (2.6)

The third term on the right-hand side describes that cells are reaching and bouncing on the vessel wall. Also,

$x_{i}^{0}$ is an initial condition of x_i.

2.1.2. Polarity of cell

We present a governing equation of a polarity q_i with the index i at a time t, defined as

$\frac{d q_{i}}{d t} = M_{q} (q_{0}^{2} - | q_{i} |^{2}) q_{i} + a_{3} \nabla c + ζ_{i} + a_{4} \frac{d x_{i}}{d t} i n Ω \times (0, T),$ (2.7)

$q_{i} = q_{i}^{0} a t t = 0.$ (2.8)

The first term on the right-hand-side has an effect of adjustment by which polarity q_i comes to meet

$q_{0}^{2} = | q_{i} |^{2}$ , where q₀ is decided arbitrarily. Chemotaxis is modeled by

$\nabla c$ and

$ζ_{i}$ represents a random number obtained under a stochastic distribution. Fictitious force is given as

$\frac{d x_{i}}{d t}$ . Also,

$q_{i}^{0}$ is an initial condition of q_i.

2.1.3. Field of fibronectin concentration

Finally, for real-valued parameters D, M, μ, a governing equation for a fibronectin concentration is defined as

$\frac{d f}{d t} = [\nabla \cdot M \nabla {D Δ f + w^{'}} - f μ] H in Ω \times (0, T),$ (2.9)

$n \cdot \nabla f = 0 o n \partial Ω \times (0, T),$ (2.10)

$n \cdot \nabla (D Δ f + w^{'}) = 0 o n \partial Ω \times (0, T),$ (2.11)

$f = f^{0} a t t = 0,$ (2.12)

where f⁰ is an initial condition for f, and

$H (x, x_{i}, t) = \sum_{i = 1}^{N_{c e l l}} e x p (- \frac{| x - x_{i} |}{a_{5}}),$ (2.13)

$w^{'} (f) = 2 f (f - 1) (2 f - 1) .$ (2.14)

3. Numerical scheme

Letting t, Δt respectively stand for time and a time increment, then from discretizing the governing equation in the time direction at a time step 0 < n ≤ N for t = nΔt, one can deduce

$\frac{x_{i}^{n + 1} - x_{i}^{n}}{Δ t} = ϕ (f^{n}) {{\hat{q}}_{i}^{n} + \sum_{j \neq i, j = 1}^{N_{c e l l}} g_{1} \frac{x_{j}^{n} - x_{i}^{n}}{r_{i, j}^{n}}} - a_{2} \nabla f^{n},$ (3.1)

$\frac{q_{i}^{n + 1} - q_{i}^{n}}{Δ t} = M_{q} (q_{0}^{2} - | q_{i}^{n} |^{2}) q_{i}^{n} + \nabla c + ζ_{i}^{n} + a_{3} \frac{x_{i}^{n + 1} - x_{i}^{n}}{Δ t},$ (3.2)

$H^{n + 1} = \sum_{i = 1}^{N_{c e l l}} e x p (- \frac{| x - x_{i}^{n + 1} |}{a_{5}}),$ (3.3)

$\frac{f^{n + 1} - f^{n}}{Δ t} = \nabla \cdot [M \nabla {D Δ f^{n + 1} + w^{'} (f^{n})} - f^{n + 1} μ] H^{n + 1} .$ (3.4)

A random number

$ζ_{i}^{n}$ is defined as

$ζ_{i}^{n} = ζ_{i}^{n - 1} - α ζ_{i}^{n - 1} + N \sqrt{Δ t},$ (3.5)

and

$N \in ℝ^{d}$ is a vector-valued function. A standard distribution with mean 0 and variance 1 is set at every element. All numerical simulations are performed by Freefem++ [26]. Especially, a fibronectin concentration is governed by the fourth-order partial differential equation. This paper using Morley element implemented in Freefem++ [26] to calculation the weak form of the Cahn–Hilliard equation.

4. Numerical result

A domain Ω is defined as

$Ω = {x = {[x, y]}^{T}; 0 \leq x \leq 50, - 25 \leq y \leq 25} .$ (4.1)

4.1. Basic features of the model

Before applying and assessing our suggested mathematical model, we confirm basic features of the model for particle collision and chemotaxis. For all cases, 10 particles are prepared initially. The position vector x_i and polarity q_i of each cell are expressed as

$x_{i} = [25 + N (0,1), N (0,1)],$ (4.2)

$q_{i} = [N (0,1), N (0,1)],$ (4.3)

where all the parameters are listed in Table 1.

Table 1. Parameters.

Symbol	Parameter	Value
a₁	Cell mobility	10
a₂	Repulsive force between cells and fibronectin	1
a₃	Chemotaxis to VEGF concentration	1
a₄	Cell alignment effect	0.1
a₅	Radius of influence of cells	0.1
D	Square of transition area width of fibronectin phase variable	1
M	Mobility of fibronectin regions	1
μ	Fibronectin Enzyme	1,10,100
q₀	Reference size of cell polarity	10
α	Persistence of cell polarity	0.1
d_core	Diameter of the cell (repulsive area)	12
d_neutral	Diameter of the cell (neutral area)	1.5 d_core
d_adh	Diameter of the cell (adhesive area)	1.534 d_core
d_reach	Diameter of the cell (interactable area)	1.834 d_core
F_rep	Intercellular repulsive force	5
F_att	Intercellular attraction force	1
Δt	Time increment	0.002

| Show Table

DownLoad: CSV

This paper is our first trial to simulate a dynamics of angiogenesis by coupling multi-agent model and phase-field model. Therefore, there is few of previous reserches, and we set most of arbitrary parameters in temporal and spatial scales. As a result, this numerical results might be helpful to understand a dynamics of angiogenesis qualitatively, not quantitatively.

4.1.1. Particle collision

For investigation of particle collisions, VEGF is ignored with c = 0. The interparticle distance r_i_,j is much less than d_core = 12 because of an initial position vector x_i = [25 + N(0,1), N(0,1)]. Therefore all particles are mutually repulsive, with no effects of chemotaxis by c = 0. Figures 1 and 2 show f and H to visualize the time evolution of vessel development and cell behaviors. In the early stages of numerical calculations, effects of repulsive force are the spreading of all cells. Subsequently, such cells with no effect of chemotaxis move slowly and become distant in the radial direction.

Figure 1. Time evolution of f without chemotaxis for c = 0.

DownLoad: Full-Size Img PowerPoint

Figure 2. Time evolution of H without chemotaxis for c = 0.

DownLoad: Full-Size Img PowerPoint

4.1.2. Chemotaxis

Next, chemotaxis generated by VEGF is considered numerically with c = 100x. In the early stages of numerical calculations, effects of repulsive force are spreading of all cells, as discussed in the section above, because of stronger repulsive force than chemotaxis. Also, interparticle distance r_i,j increases. After the repulsive force weakens sufficiently to become weaker than chemotaxis, all particles are moving in the right direction, where VEGF c has a larger value, as portrayed in Figures 3 and 4.

Figure 3. Time evolution of f with chemotaxis for c = 100x.

DownLoad: Full-Size Img PowerPoint

Figure 4. Time evolution of H with chemotaxis for c = 100x.

DownLoad: Full-Size Img PowerPoint

4.2. Numerical simulations of angiogenesis

Based on our suggested mathematical model and numerical techniques constructed in the sections presented above, we are performing numerical simulations not only of lateral inhibitation but also of elongation, branching, which are basic features of angiogenesis. The initial condition of position vector x_i and polarity q_i of each cell is expressed as

$x_{i} = [0.05 + N (0,1), N (0,1)],$ (4.4)

$q_{i} = [1,0],$ (4.5)

where the number of cells is 20. Then we set x_i = 0.05 in the case of x_i < 0; also, c = 100x is used as VEGF. In fact, as a governing equation of fibronectin concentration f, the Cahn–Hilliard equation is representative of PDE models, but reduces to an ODE model by substituting M = 0. We calculate both to compare the results obtained from calculations.

Figure 5 portrays the time evolution of fibronectin f obtained in the case of ODE model with M = 0. The model has behavior of lateral inhibition and elongation, branching, and anastomosis. By contrast, Figure 6 portrays results obtained for the PDE model. The Cahn–Hilliard equation shows fibronectin f to be a distribution from 0 to 1. The structure of the vessel is therefore visualized as clearer than that of the ODE model. Next, we present Figure 7 in which the x coordinate of all cells is shown. As Figure 7 shows, all cells are moving quickly in the early stage because of mutual repulsion. Thereafter, chemotaxis ∇c attracts them to move slowly in the positive-x direction. Figure 8 presents data for comparison, illustrating differences of the same particles between the ODE model and PDE model. Some particles are observed to move with very different behaviors because of the effects of distribution of fibronectin f attributable to the distribution of the fibronectin concentration f in the ODE model, which differs from that of the PDE model.

Figure 5. Time evolution of f with chemotaxis in the case of c = 100x for the ODE model.

DownLoad: Full-Size Img PowerPoint

Figure 6. Time evolution of f with chemotaxis for c = 100x for the PDE model.

DownLoad: Full-Size Img PowerPoint

Figure 7. x coordinates of all cells with time steps up to t = 20.

DownLoad: Full-Size Img PowerPoint

Figure 8. Particle difference between M = 0 and M = 1.

DownLoad: Full-Size Img PowerPoint

4.3. Attraction of vessels by tumor cells

Finally, on the assumption that tumor cells express vascular endothelial growth factor receptors and respond to autocrine and paracrine VEGF signals, this section defines VEFG as

$c = 1000 e x p [- 0.001 {{(x - 50)}^{2} + y^{2}}],$ (4.6)

where Figure 9 depicts the distribution of this c. An initial condition of position vector x_i and polarity q_i of each cell is expressed as

$x_{i} = [0.05 + N (0,1), N (0,6)],$ (4.7)

$q_{i} = [1,0],$ (4.8)

where the number of cells is 20. We set x_i = 0.05 in the case of x_i < 0. By the way, N(0,1) is used in section 4.1 and 4.2. In the case of N(0,1), it is seemed that effects of cross-interaction between particles themselves are strong in eraly stage. Thereofre, we set N(0,6) in section 4.3 to focus on effects of a tumor cell.

Figure 10 shows the time evolution of f with chemotaxis at t = 62, unlike in the case of Figures 5 and 6, not all particles are moving in the x-direction. They coming in the direction of a coordinate (x,y) = (50,0) with a maximum value of c. However, apparently, some particles are going to the lower right. Therefore we are concerned that this phenomenon is attributable to cell adhesion effects. The domain used for this simulation is 50 × 50 square. The cell diameter is d_adh ≈ 18 from Table 1. Therefore, the top position cell affected by VEGF is moving in the lower right. Subsequent cells are following the top position cell because of cell adhesion.

Next, we observe vessel development depending on μ representing the strength for degenerating fibronectin and its corresponding x coordinate at each cell in the case of the PDE model. As presented in Figures 11–13, a small μ as in Figure 11 with μ = 1 can suppress vessel development more than that shown in Figure 13 with μ = 100 because a cell can be moving much more in the area of fibronectin f = 0 than f = 1. Consequently, numerical simulations of Figures 11–13 represent biochemically adequate numerical results.

Figure 9. Distribution of c.

DownLoad: Full-Size Img PowerPoint

Figure 10. Time evolution of f.

DownLoad: Full-Size Img PowerPoint

Figure 11. (a) Fibronectin f at t = 30 and (b) x coordinate at each cells up to t = 20 for μ = 1.

DownLoad: Full-Size Img PowerPoint

Figure 12. (a) Fibronectin f at t = 30 and (b) x coordinate at each cells up to t = 20 for μ = 10.

DownLoad: Full-Size Img PowerPoint

Figure 13. (a) Fibronectin f at t = 30 and (b) x coordinate at each cells up to t = 20 for μ = 100.

DownLoad: Full-Size Img PowerPoint

5. Conclusion

We demonstrated the characteristic features of angiogenic endothelial cell (EC) behaviors. Taken together with findings reported from earlier studies using experimentation such as that reported by Arima et al. [18], we suggest a multi-agent model coupled with phase-field model for angiogenic morphogenesis. Such a model enables us to dissect cellular mechanisms regulating dynamic and complex multicellular processes systematically. The Cahn–Hilliard equation is regarded as phase-field model, and ODE and PDE equation are prepared by switching M = 0,1.

Our first numerical result is that multi-agent based EC movements with ODE model, represented as a process driven by a simple stochastic rule in this suggested multi-agent model, can adequately describe branch elongation, which plays an important role as one angiogenic model. However, EC behaviors were insufficiently explained using cell processes alone because the ODE model was not able to consider an effect of elasticity of a vessel wall. This numerical result leads to our second result suggesting that PDE model might be necessary for degenerating fibronectin f in the whole cell. This prediction was verified biologically. Consequently, this study advances our mathematical understanding of angiogenic morphogenesis as coordinated multicellular processes.

As described herein, a smaller area is adopted to suppress calculation cost to the greatest degree possible. Effects of cell adhesion are increased as result. It is too difficult to survey some biochemical phenomena precisely. Future work must be done for numerical calculation with a larger area than this study to investigate additional details about relations with chemotaxis ∇c and enzyme of fibronectin μ.

In fact, we need to analyze correlation between chemotaxis and heterogeneous directionality, cell mixing, cell overtaking, branching morphogenesis. However, it is seemed that effects of chemotaxis is much stronger than that of them, and we did not observe about details of the movement of ECs from our numerical simulation.

In the future, a multi cellular model suggested here will be quantitatively leading a life science reserch on molecular level and a cellular level to acquisitions of a macroscopic histogenesis and a vital function.

Acknowledgments

This work was supported by JSPS KAKENHI, Grant Number 19K03645, and MEXT KAKENHI, Grant Number 17H06327.

Conflict of interest

The authors declare no conflicts of interest in this paper.

References

[1]	Holderfield MT, Hughes CC (2008) Crosstalk between vascular endothelial growth factor, notch, and transforming growth factor-beta in vascular morphogenesis. Circ Res 102: 637-652. https://doi.org/10.1161/CIRCRESAHA.107.167171
[2]	De Smet F, Segura I, De Bock K, et al. (2009) Mechanics of vessels branching: filopodia on endothelial tip cells lead the way. Arterioscle Thromb Vasc Biol 29: 639-649. https://doi.org/10.1161/ATVBAHA.109.185165
[3]	Armulik A, Abramsson A, Betsholtz C (2005) Endothelial/pericyte interactions. Circ Res 97: 512-523. https://doi.org/10.1161/01.RES.0000182903.16652.d7
[4]	Gaengel K, Genove G, Armulik A, et al. (2009) Endothelial-mural cell signaling in vascular development and angiogenesis. Arterioscle Thromb Vasc Biol 29: 630-638. https://doi.org/10.1161/ATVBAHA.107.161521
[5]	Hellstrom M, Gerhardt H, Kalen M, et al. (2001) Lack of pericytes engenders endothelial hyperplasia and abnormal vascular morphogenesis. J Cell Biol 153: 543-554. https://doi.org/10.1083/jcb.153.3.543
[6]	Lafleur MA, Forsyth PA, Atkinson SJ, et al. (2001) Perivascular cells regulate endothelial membrane type-1 matrix metalloprotainase activity. Biochem Biophy Res 282: 463-473. https://doi.org/10.1006/bbrc.2001.4596
[7]	Liu H, Kennard S, Lilly B (2009) NOTCH3 expression is induced in mural cells through an autoregulatory loop that requires endothelial-expressed JAGGED1. Circ Res 104: 466-475. https://doi.org/10.1161/CIRCRESAHA.108.184846
[8]	Frieboes HB, Lowengrub JS, Wise SM, et al. (2007) Computer simulation of glioma and morphology. Neuroimage 37: 59-70. https://doi.org/10.1016/j.neuroimage.2007.03.008
[9]	Wise SM, Lowengrub JS, Frieboes HB, et al. (2008) Three-dimensional multispecies nonlinear tumor growth – 1: Model and numerical method. J Theor Biol 253: 524-543. https://doi.org/10.1016/j.jtbi.2008.03.027
[10]	Frieboes HB, Jin F, Chuang YL, et al. (2010) Three-dimensional multispecies nonlinear tumor growth – 2: tumor invasion and angiogenesis. Neuroimage 264: 1254-1278. https://doi.org/10.1016/j.jtbi.2010.02.036
[11]	Lima EA, Oden JT, Almeida RC (2014) A hybrid ten-species phase-field model of tumor growth. Math Models Methods Appl Sci 24: 2569-2599. https://doi.org/10.1142/S0218202514500304
[12]	Perfahl H, Bryne HM, Chen T, et al. (2011) Multiscale modelling of vascular tumor growth in 3D: The role of domain size and boundary conditions. PloS One 6: e14790. https://doi.org/10.1371/journal.pone.0014790
[13]	Alber MS, Kiskowski MA, Glazierz JA, et al. (2003) On cellular automaton approaches to modelling biological cells. Mathematical Systems Theory in Biology, Communications, Computation, and Finance. New York: Springer 1-39. https://doi.org/10.1007/978-0-387-21696-6_1
[14]	McAneney H, O'Rourke SFC (2007) Investigation of various growth mechanisms of solid tumor growth within the linear-quadratic model for radiotherapy. Phys Med Biol 52: 1039-1054. https://doi.org/10.1088/0031-9155/52/4/012
[15]	Rejniak KA, Anderson ARA (2010) Hybrid models of tumor growth. Wires SystBiol Med 3: 115-125. https://doi.org/10.1002/wsbm.102
[16]	Boondirek A, Lenbury Y (2006) A stochastic model of cancer growth with immune response. J Kor Phy Phys Soc 4: 1652-1666.
[17]	Gerlee P, Anderson ARA (2008) A hybrid cellular automaton model of clonal evolution in cancer: The emergence of the glycolytic phenotype. J Thero Biol 250: 705-722. https://doi.org/10.1016/j.jtbi.2007.10.038
[18]	Arima S, Nishiyama K, Ko T, et al. (2011) Angeogenesis morphogenesis driven by dynamic and heterogeneous collective endothelial cell movement. Development 138: 4763-4776. https://doi.org/10.1242/dev.068023
[19]	Takubo N, Yura F, Naemura K, et al. (2019) Cohesive and anisotropic vascular endotherial cell motility driving angiogenic morphogenesis. Scientific Reports 9: 9034. https://doi.org/10.1038/s41598-019-45666-2
[20]	Cahn JW, Hilliard JE (1958) Free energy of a nonuniform system. I. Interfacial free generic interface. J Chem Phys 28: 258-267. https://doi.org/10.1063/1.1744102
[21]	Furihata D, Matsuo T (2011) Discrete Variational Derivative Method: A structure-preserving numerical method for partial differential equations. Nwe York: CRC Press. https://doi.org/10.1201/b10387
[22]	Feng X, Li Y, Xing Y (2016) Analysis of mixed interior penalty discontinuous Galerkin methods for the Cahn–Hilliard equation and the Hele–Shaw flow. SIAM J Numer Anal 54: 825-847. https://doi.org/10.1137/15M1009962
[23]	Kovacs M, Larsson S, Mesforsh A (2011) Finite element approximation of the Cahn–Hilliard–Cook equation. SIAM J Numer Anal 49: 2407-2429. https://doi.org/10.1137/110828150
[24]	Brenner SC, Sung Li-y, Zhang H, et al. (2013) A Morley finite element method for the displacement obstacle problem of clamped Kirchhoff plate. J Comp Appl Math 254: 31-42. https://doi.org/10.1016/j.cam.2013.02.028
[25]	Li M, Guan X, Mao S (2014) New error estimates of the Morley element for the plate bending problems. J Comp Appl Math 263: 405-416. https://doi.org/10.1016/j.cam.2013.12.024
[26]	Hecht F (2012) New development in FreeFem++. J Numerical Math 20: 251-266. https://doi.org/10.1515/jnum-2012-0013
[27]	Sánchez-Pérez R, Pavan S, Mazzeo R, et al. (2018) Mutation of a bHLH transcription factor allowed almond domestication. Science 364: 1095-1098. https://doi.org/10.1126/science.aav8197

This article has been cited by:

1.	Egor V. Yakovlev, Ivan V. Simkin, Anastasiya A. Shirokova, Nataliya A. Kolotieva, Svetlana V. Novikova, Artur D. Nasyrov, Ilya R. Denisenko, Konstantin D. Gursky, Ivan N. Shishkov, Diana E. Narzaeva, Alla B. Salmina, Stanislav O. Yurchenko, Nikita P. Kryuchkov, Machine learning approach for recognition and morphological analysis of isolated astrocytes in phase contrast microscopy, 2024, 14, 2045-2322, 10.1038/s41598-024-59773-2
2.	Veronika Kopylova, Stanislav Boronovskiy, Yaroslav Nartsissov, Approaches to vascular network, blood flow, and metabolite distribution modeling in brain tissue, 2023, 15, 1867-2450, 1335, 10.1007/s12551-023-01106-0

Reader Comments

Your name:*

Email:*
© 2022 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)

通讯作者: 陈斌, bchen63@163.com

1.
沈阳化工大学材料科学与工程学院沈阳 110142

AIMS Bioengineering

Metrics

Article views(2315) PDF downloads(156) Cited by(2)

Preview PDF

Download XML

Export Citation

Article outline

Show full outline

Figures and Tables

Figures(13) / Tables(1)

AIMS Bioengineering

Multicellular model of angiogenesis

Related Papers:

Abstract

1. Introduction

2. Problem formulation

2.1. Governing equation

2.1.1. Equation of cell motion

2.1.2. Polarity of cell

2.1.3. Field of fibronectin concentration

3. Numerical scheme

4. Numerical result

4.1. Basic features of the model

4.1.1. Particle collision

4.1.2. Chemotaxis

4.2. Numerical simulations of angiogenesis

4.3. Attraction of vessels by tumor cells

5. Conclusion

References

This article has been cited by:

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Figures and Tables

Other Articles By Authors

Catalog

AIMS Bioengineering

Multicellular model of angiogenesis

Related Papers:

Abstract

1. Introduction

2. Problem formulation

2.1. Governing equation

2.1.1. Equation of cell motion

2.1.2. Polarity of cell

2.1.3. Field of fibronectin concentration

3. Numerical scheme

4. Numerical result

4.1. Basic features of the model

4.1.1. Particle collision

4.1.2. Chemotaxis

4.2. Numerical simulations of angiogenesis

4.3. Attraction of vessels by tumor cells

5. Conclusion

References

This article has been cited by:

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Figures and Tables

Other Articles By Authors

Related pages

Tools

Export File

Citation

Format

Content

Catalog