Modelling of the migration of endothelial cells on bioactive micropatterned polymers

: In this paper a macroscopic model describing endothelial cells migration on bioactive micropatterned polymers is presented. It is based on a system of partial diﬀerential equations of Patlak-Keller-Segel type that describes the evolution of the cell densities. The model is studied mathematically and numerically. We prove existence and uniqueness results of the solution to the diﬀerential system and also that fondamental physical properties such as mass conservation, positivity and boundedness of the solution are satisﬁed. The numerical study allows us to show that the model behaves in good agreement with the experiments. Résumé : Dans cet article nous présentons un modèle macroscopique décrivant la migration de cellules endothéliales sur un micropattern de polymers bioactifs. Ce modèle est basé sur un système d’équations aux dérivées partielles du type Patlak-Keller-Segel. Les propriétés mathé-matiques et numériques du modèles sont présentées. Nous démontrons des résultats d’existence et d’uncité ainsi que les propriétés physiques, telles que la conservation de la masse, la positivité et le caractère borné de la solution. L’étude numérique nous permet de montrer que le modèle corrobore les résultats expérimentaux.


Introduction
Tissue engineering is the use of combination of cells, engineering, materials, and suitable biochemical factors to improve or replace biological functions [26]. Tissue engineering has being quickly developing since the 1990s [26]. However, the major roadblock for engineering large tissues is the lack of functional microvasculature networks, which provide nutrients and oxygen for tissue survival and remove the waste product from metabolism [18]. The lack of vascularization has hampered the survival of engineered tissues after implantation [18]. Researchers rely on the increasing knowledge of angiogenic and vasculogenic processes to stimulate vascular network formation [32,31]. This complex process of new blood vessel formation is orchestrated by the interaction between endothelial cells (ECs) and their neighboring mural cells via a complex network of intracellular signaling mechanisms [28,17]. Ever since the introduction of the in vitro model of angiogenesis [11], there has been great research interest to understand the intricate process of tube formation. Although many efforts have been made, the mechanism associated with angiogenesis and vascularization is poorly understood. A deeper comprehension of cellsbiomaterials interaction is required for basic understanding of angiogenesis and vascularization in tissue engineering [5].
One strategy in developing clinical implants is therefore to use bioactive materials which can either elicit a regenerative response at the site of damage in vivo or be used to grow tissue in vitro for subsequent implantation [2,23]. Different bioactive ligands have been used to study their effects on cell functions for a better understanding of vascularization [31]. In the aim at promoting angiogenesis in the case of tissue engineering or inhibiting angiogenesis in the case of cancer, it is now important to understand the mechanisms that regulate lumen formation. Successful micropatterning of cells is becoming a key component of this field [16]. Researchers are now interested in the behavior of cells on substrates that have been patterned by micro/ nanofabrication [10,27]. It is known that cell positioning and physiology can be controlled by the substrate on which the cells adhere [6]. In our study, using cell adhesion peptides micropatterned onto material, we observed a tube-like formation in comparison to virgin material or homogeneously grafted materials [22,23].
Experimental studies using micropatterned substrates revealed that the cell migration is governed by the geometry of the patterns (size of patterns and distance between patterns). Endothelial cells so cultured form extensive cell-cell interactions. Accumulation of endothelial cells junctions implies that some cells form tube-like structures. The goal of the present paper is to provide a model that describes such results.
Adhesive areas are composed of cell adhesion peptides or growth factors peptides that make the cells adhere. These areas are surrounded by non-adhesive areas [22]. We assume (and this is actually confirmed by experiments) that active principles (cell adhesion peptides or growth factors) are not diffusing spatially. Therefore, endothelial cells outside adhesive areas have no mean to "feel" directly these areas. They find out these areas indirectly. We do not consider the influence of nutrients and assume that cells obtain enough nutrients from the material (due to grafted active principles onto material). Endothelial cells are seeded onto micropatterned bioactive materials during several hours, then they are washed out. Only the adhered endothelial cells remain on material. The initial density of cells is around 40 000 cells per cm 2 . At the beginning of the experiments, during the migration phase we observe that cells have a random motility and stop on adhesive areas. Moreover the attraction of endothelial cells on adhesive areas seems to be higher than the one of cells outside these areas. Experiments show that endothelial cells are grouping together along the micropatterns. According to pattern size, endothelial cells line their cytoskeleton to adjust it with the adhesive area. One can also notice that, with micropatterns of 10 µm thin strips, tubes containing a central lumen may appear [23,7]. In Inria other words, blood vessels are created from an initial random density of endothelial cells. Such phenomenon is not observed with larger strips [7,20,25].
To illustrate these experiments, we present in Fig.1 pictures of the micropatterned bioactive materials at the end of the migration step. Two different micropatterns are considered: on Fig.1(a) thin adhesive areas (bioactive pattern size: 10 µm and distance between patterns: 100 µm) have beed used, whereas Fig.1(b) shows the end of the migration on large strips (bioactive pattern size: 300 µm and distance between patterns: 100 µm). For a detailed description of such experiences, one may refer to [22,23].   [22]. The distance between bioactive patterns is 100 µm.
We observed that for the largest adhesive areas the adhered cells density is smaller than for thin strips. Therefore, the geometries of the micropatterns play a crucial role in the endothelial cells migration and then in the formation of new vessels.
In this paper, we are interested in understanding how these patterns (size and distance between microfeatures) do influence endothelial cells migration. The model we present here is a continuous Patlak-Keller-Segel type model [1,13,21,30]. The chemotaxis term takes the cellcell interactions into account instead of the cell-chemical attractant interactions. We show that this new model based on a system of coupled partial differential equations conserves the mass and existence and uniqueness results of weak solution hold. We also provide numerical results in accordance with the experiments, which ensures the validity of our model. Moreover, these simulations allow us to obtain informations on the influence of the geometry and of the initial concentration of cells on the migration of the cells.
The outline of the article is the following. In section 2, we describe the mathematical model and we state the main result of global existence and uniqueness of the weak solution to the P.D.E system. Section 3 is devoted to the proof of the main theorem. We then provide numerical results in section 4 in order to compare the simulations to the experiments.

Description of the model and main result
In this section, we describe the Patlak-Keller-Segel type continuous model we study throughout the paper. The model is composed of a diffusion term coupled with a reaction term, which describes the effect of the chemoattractants, which themselves statisfy a diffusion equation.
Various continuous models of Patlak and Keller-Segel type have been used to describe cell motility. Typically, the governing equations of these models are written in the following general form: on a domain Ω ⊂ R n where u denotes the cells density, v is the chemical signal concentration. The diffusive terms take the random motility of the cells into account, whereas the advection describes the influence of the chemical signal on the motion of the cells. The two corresponding diffusion parameters are denoted by D 1 and D 2 while χ is the chemotaxis coefficient. The function f describes the growth and the death of cells, whereas g and h are the production and degradation of the chemotaxic signal. These equations have been theoretically studied [3,4,8,12,14,33] for several years. Based on this extensive literature we provide a slightly modified model to describe the cell migration on bioactive micropatterns.

Statement of the equations
According to the experiments, the behavior of the cells is drastically different on the adhesive area and outside this area. Actually, outside the adhesive strip the cells seem to attract each other (probably thanks to the chemoattractant they produce) and also diffuse in the domain, but as soon as they reach the adhesive strip they seem stuck on the strip and then they diffuse only on this area, ignoring the outer cells. Moreover it seems that the production of chemoattractant of the cells located on the adhesive strip is bigger than outside.
Since there is no clear understanding of the way that endothelial cells communicate, we chose to consider the chemotaxis term as the attraction between endothelial cells (and we do not consider any gradient of concentration of the chemoattractant).
Based on these assumptions, we derive the following model. Consider a domain Ω composed by the adhesive areas, denoted Ω, and the non-adhesive areas Ω\ Ω, both domains being bounded domains with smooth boundary.
Two different types of endothelial cells are considered. We denote by u 1 (t, x, y) the density of endothelial cells, at any point (x, y) and time t, that can freely move (i.e. they have yet to move over adhesion proteins). Cells that are adhering on the substrate are tracked through their density u 2 . The function v represents the density of the chemoattractant. The equations governing the endothelial cells migration are given for t > 0 by in Ω, (1a) with the homogeneous boundary conditions on ∂Ω and ∂ Ω: and with the initial conditions (u 0 1 , u 0 2 , 0) We then denote by u the sum The parameters d 1 , d 2 , η, γ 1 , γ 2 and λ are strictly positive and they will be fitted by the experiments in a forthcoming work, but we consider here that they are given constants. The coefficients d 1 and d 2 denote the diffusion coefficients of the cells u 1 and u 2 respectively. The coefficient η > 0 is the self-degradation rate of the chemoattrant produced by the cells, while the coefficients γ 1 and γ 2 are the coefficients of the production of the chemoattractant respectively for the cell u 1 and u 2 . The parameter λ is the rate for the cell u 1 to become u 2 as soon as u 1 lies in the micropatterns. The two first equations describe the migration of the cells on Ω. Outside the strips, endothelial cells diffuse and attract the neighboring cells via the chemotaxis sensitivity function : Here above, χ 0 is a chemotaxis parameter, and the term (1 − u 1 ) is settled to prevent the overcrowding of the cells u 1 . Endothelial cells once they reach the adhesive area Ω are captured and then diffuse only in the strip. This is handled by the penalty term −λ1 Ω u 1 (1 − u 2 ). Cells on the strips still have a random motility and their concentration grows up as the term λ1 Ω u 1 (1−u 2 ), where 1 − u 2 prevents the blow-up of u 2 in equation (1b). The third equation describes the production of the chemoattractant by the cells. Since the cells on the strip seem to be more attractive, we suppose that the production coefficients checks 0 < γ 1 < γ 2 . We also add a degradation coefficient η > 0 describing the metabolization of the chemoattractant.

Main theoretical result
We have the following theorem which is a straightforward consequence of the results of the next section 3: Theorem 2.1. Let d 1 , d 2 , η, γ 1 , γ 2 and λ be strictly positive constants. Suppose that the initial data There exists a unique weak solution (u 1 , u 2 , v) to problem (1) such that and for almost any t > 0 The next section is devoted to prove this theorem. The proof is based on Gaussian upper bounds for heat kernels [29]- [35].

Theoretical study of the model
In this section, we study the mathematical properties of the model. Throughout this section we suppose that Ω and Ω are smooth domains of R 2 . We remind that d 1 , d 2 and η are strictly positive coefficients.

Kernels of the operators
The aim of this paragraph is to define and provide estimates of the kernels of the operators ∂ t − ∆ + η and ∂ t − d 1 ∆ on Ω and for the kernel of ∂ t − d 2 ∆ on Ω, with homogeneous Neumann conditions respectively on ∂Ω and ∂ Ω.
on Ω, all with homogeneous Neumann conditions, are respectively defined by and for any (t, y) ∈ (0, ∞) × Ω, for B, while G is given by and for any (t, y) ∈ (0, ∞) × Ω, and G is the solution to and for any (t, y) ∈ (0, ∞) × Ω, Note that the above kernels are symmetric in their second and third variables.
Proposition 1. For any y ∈ Ω (respectively for any y ∈ Ω), we have the following estimates for positive constants C Ω and C Ω , which depend on the domain Ω and Ω respectively: and gradient estimates hold too:

Inria
In addition, due to the boundedness of Ω, we also have ∇ y G(t, ·, y) Proof. Obviously the coefficients of diffusion d 1 and d 2 , since they are strictly positive constants do not play a crucial role, and can be supposed to be 1, after an appropriate rescaling of the time variable t. Moreover it is sufficient to prove the above estimates for the heat kernel G, since For t ≥ 1, estimates (3.2)-(3.3) of [35] straightforwardly provide the result. Suppose that 0 < t ≤ 1. Estimates (5a) easily come from Theorem 6.10 pp 171 of [29], since for any x ∈ Ω, Estimates (6a) are consequences of the section 6.6 entitled Weighted Gradient Estimates and in particular Theorem 6.19 p 185 [29]. Actually by Cauchy-Schwarz inequality hence the estimates (6a). Now let φ ∈ L ∞ (Ω), by estimates (6a) and since the measure |Ω| of Ω is bounded we infer Ω |φ(y)| Ω |∇ y G(t, x, y)| dx dy = Ω×Ω |φ(y)| |∇ y G(t, x, y)| dy dx, hence estimates (7a), which ends the proof of the proposition.
Remark 1. The above estimates are probably not optimal, since for the half-plane the heat kernel write: and therefore the power t −3/4 has to be replaced by t −1/2 similarly to the heat kernel of the whole plane R 2 . However these results are sufficient to prove existence and uniqueness of the solution to problem (1).

Inria
Using estimates (5a)-(6a) we infer that for any (ν 1 , ν 2 ) ∈ X T M × X T M : Define now the operator T on X T M × X T M by where T 1 is the operator defined on X M × X T M by and T 2 is defined by Remark 2. Proving that T is a contraction mapping from X T M ×X T M onto itself for small enough time T will then ensure the local existence of the weak solution given by (10) to problem (1).

Contraction mappings
Proposition 2. The operators T is a contraction mapping from X T M × X T M onto itself for T small enough.
Proof. The proof is based on the properties of the kernel B, G and G given by Proposition 1. Thanks to estimate (5a) we deduce for any (ν 1 , ν 2 ) ∈ X T M × X T M : hence for T small enough T 2 maps X T M × X T M onto X T M . Moreover using inequality: we infer for T small enough the operator T 2 is a contraction mapping from X T M × X T M onto X T M . Prove now that T 1 is a contraction mapping from X T M × X T M onto X T M . First observe that for any s ∈ R, |s|/(1 + |s|) ≤ 1 hence for any ν 1 ∈ X T M , for any s ∈ R, hence for any (ν 1 , ν 2 ) ∈ X T M × X T M χ 1 (ν 1 , L(ν 1 , ν 2 ))(t, ·) L ∞ (Ω) ≤ χ 0 (1 + M ), for almost any t ∈ (0, T ), and thanks to estimates (6a)-(11) This implies that for T small enough T 1 maps X T M × X T M onto X T M . In addition, observe that for two couples (ν 1 , ν 2 ) and (µ 1 , µ 2 ) belonging to X T M × X T M we have where to simplify notations we have denoted by χ ν1,ν2 the function and similarly for χ µ1,µ2 . According to estimate (6a), and thanks to the definition of L, we infer Moreover, observing that we deduce from estimates (5a)-(6a)-(7a) and (11) that there exists a constant C > 0 which depends on M , and on the parameters χ 0 , γ 1 , γ 2 , λ such that which ensures the strict contractility of T 1 for T small enough, and therefore T is a strict contraction from X T M × X T M onto itself.
The Picard fixed point theorem straightforwardly implies the following theorem of existence and uniqueness for small time.
Then, for T small enough, there exists a unique weak solution (u 1 , u 2 , v) to (1) such that

Mass conservation and global existence
We first observe that the total mass of cells is conserved.
and let T small enough so that a weak solution (10) to (1) exists. Then, for any t ∈ [0, T ] we infer Proof. Actually integrating (1a) and (1b) respectively and summing the integrands imply, since ∂ n u 1 | ∂Ω , ∂ n u 2 | ∂ Ω and ∂ n v| ∂Ω vanish We now show that if u 0 1 and u 0 2 are positive and bounded by 1 then u 1 and u 2 stay positive and bounded by 1 on [0, T ].
and let T small enough so that a weak solution given by (10) to problem (1) exists. If (u 0 1 , u 0 2 ) are such that In addition, Therefore, the weak solution (10) exists for almost any t ∈ (0, +∞).
Proof. First observe that if u 1 is positive then since u 0 2 is positive the function u 2 is positive almost everywhere. Actually multiplying (1b) by u − 2 = max(0, −u 2 ) and integrating by parts implies hence u − 2 equal zero by Gronwall's lemma. Similarly, if u 1 is positive, then since u 2 is therefore also positive we infer that v is positive by multiplying (1c) by v − and integrating by parts.
Similarly let Then, U 1 satisfies Once again multiply (12) by U + 1 = max(U 1 , 0) and integrate by parts to obtain Since 1 − u 2 is positive and using Cauchy-Schwarz estimate and Peetre inequality for α large enough (as used above to prove that u 1 ≥ 0) implies that Therefore, Gronwall lemma implies that U + 1 vanishes almost everywhere in (0, T ) × Ω hence u 1 ≤ 1.
To obtain the positivity of v, first multiply (1c) by v − and integrate by part to infer, since u 1 and u 2 are positive that: Since γ 1 (u 1 − 1) + γ 2 (u 2 − 1) ≤ 0, we infer that V + identically vanishes after multiplication and integration by parts, hence From the implicit representation integral of u 1 and u 2 we deduce easily that if T M is the maximal time of existence, then there exists a sequence (t n ) n∈N tending to T M , with t n < T M such that lim n→+∞ u 1 (t n , ·) L ∞ = +∞, hence u 1 and u 2 exists for almost any t ∈ (0, +∞) by contraposition.

Inria
Theorem 2.1 is an easy consequence of the above results.
The following result is a straightforward consequence of proposition 4. It ensures that the mass of the cells tends to concentrate on the micropatterns. and let (u 1 , u 2 ) the weak solution to problem (1). Then

Numerical results
We provide now numerical schemes that are used to compute problem (1), and then we show the simulations that corroborate the experimental results.

Approximation of the problem
We consider a cartesian mesh (composed by quadrilaterals). We discretize the model using finite volume method [9] and we use an implicit Crank-Nicolson scheme for the time discretization. We solve the model using a decoupled approach [15]. In particular, the first equation is split into advection and diffusion parts. Let us recall the expression of this equation : To simplify the notations we define A and B as: Let us denote the time step by ∆t, set t n = n∆t and let (u n 1 , u n 2 , v n ) be the solution at the time t n . At each time step we first solve the diffusive part : . For all the diffusive terms, the spatial discretization is handled by a centered finite volume scheme, all the species being computed at the centre of each element of the mesh. We then solve the advection part : The high order WENO 5 (Weighted Essentially Non-Oscillatory) finite difference scheme introduced in [24] and improved in [19] has been used to approach the convective term. These solvers are implemented in the academic library eLYSe 1 . In the following, the initial conditions write where u 0 is a function of x ∈ Ω. Hence at the initial time the support of u 1 and u 2 are disjoint.

Mathematical behavior of the model
In this paragraph, we present numerical results illustrating the good mathematical behavior of our model. We want to check the properties of the model, when the maximal density on the adhesive area is reached. We consider a domain Ω = [0, 1] × [0, 1] and meshes composed by 100 × 100 quadrilaterals. The domain is composed by a unique adhesive area (cf Fig. 2).  Fig. 4. When considering u 0 = 0.08, the maximal density on the adhesive area is never reached. We observe that u 1 is decreasing, while u 2 is increasing with respect to the time. In the second case, u 0 = 0.25, the maximal density for u 2 is reached for t = 0.3. As expected the migration stops and no more movement are observed. These simulations show that a minimum amount of endothelial cells is needed at initial time in order to reach an optimal concentration on the strips at the end of the experiment. If this initial concentration is too small the final density of endothelial cells is suboptimal.

Behavior on realistic benchmarks
We now provide simulations in realistic confirgurations. Therefore, throughout this subsection the function u 0 of (14) is a random spatial distribution between 0 and 1 using a normal distribution.

Behavior on thin strips
We first consider a micropattern composed by six thin adhesive areas (in red on Fig.5(b)).
The figure Fig. 7 shows the behavior of v for the same set of parameters.    We obtain a good agreement with the expected evolution. Indeed, the density of the cells on the adhesive areas increases in time, whereas outside it becomes very small. Cells are stucked on the strips and stop moving once they are over them. As a consequence the density of the attractant on the strips also increases.
In Fig. 10, we present the behavior of v for the same choice of parameters at time t = 0.3 and t = 0.6.
As previously we observe a behavior in good agreement with the experiments. When consid-    ering two large adhesive areas the velocity of the migration is smaller than for a large number of thin strips. Indeed, at the time step t = 1.0 we observe that with thin strips the migration seems to be more advanced than in the case with large strips. This could be explained by the fact this last case some cells are far away from a strip and their migration toward the strips take more time.

Influence of the number of strips on the migration
We want to study the influence of the geometry on the migration. We set the surface of the adhesive domain, and let the number of strips, N s , vary. The average of u 2 in term of the time for N s = 1, 2, and 4 is presented in Fig. 11 . We observe that when considering four strips the migration is quicker. Moreover the mean density reached is higher, which corroborates the experiments.

Conclusion
In this paper, a macroscopic model for endothelial cells migration is presented. Its major biological assumption is that cells are not actively attracted by the adhesion patterns but just adhere on it and try to gather.
Mathematically, mass conservation and global existence is shown. Numerically, the model behaves in good agreement with the biological knowledge. Despite the lack of active attraction by the adhesion patterns, the non-washed out endothelial cells end up on the patterns. We have observed two facts that have been reported by the experiments: 1) For a given surface of active principle the process of cell migration is more efficient with a large number of thin strips than with a small number of large strips.
2) There exists a minimum value of the initial density of endothelial cells to impose in order to have an optimal cell migration towards the active principle.