Some Numerical Results on Chemotactic Phenomena in Stem Cell Therapy for Cardiac Regeneration

: Biological models for cardiac regeneration and remodeling, along with the effects of cytokines or chemokines during the therapy with mesenchymal stem cells after a myocardial infarction, are of crucial importance for understanding the complex underlying mechanisms. This paper presents a mathematical model composed of three coupled partial differential equations that describes the dynamics of stem cells, nutrients and chemokines, highlighting the fundamental role of the chemokines during the myocardial tissue regeneration process. The system is solved numerically using mimetic difference operators and the MOLE library for MATLAB. The results show the tissue regeneration process in the necrotic part closest to the cell implantation area.


Introduction
In the field of regenerative medicine, stem cell therapy has been used to treat heart diseases after a myocardial infarction [1,2].Different stem cell populations have been explored for their ability to renew themselves and differentiate into multiple lineages [3,4].These treatments take into account different protocols.In particular, the number of cells implanted, the route of administration and the time of application are variable and potentially very important factors in determining the effect of the treatment [5][6][7].
Developing a mathematical model for cardiac tissue regeneration therapy with stem cells is a challenging task due to the intricate interconnections among various complex processes.In the literature, numerous studies propose mathematical models that describe only some of the processes involved, not only in myocardial-regeneration-related works but also in studies concerning the use of stem cells in tumor treatment or wound healing [8].
In this paper, a mathematical model based on a system of evolutionary and diffusive partial differential equations is proposed.The unknowns of this system will be naturally the concentration of stem cells but also those of nutrients (namely, oxygen) and messengers (chemokines).This mathematical scheme explicitly takes into account the three-dimensional geometry of the tissues, where the injured area is represented as a portion of the spatial domain in which the problem is posed.The boundary separating this area from healthy tissue moves during reconstruction and, of course, the objective is to advance toward the necrotic area as much as possible.Furthermore, the explicit spatial dependence of our scheme allows us to directly and realistically introduce, for example, the different routes of administration and quantity of implanted stem cells, which can be described by different boundary conditions.
To achieve the description of the dynamics of stem cells and the contraction of the necrotic core, we employ diffusion-reaction equations, which are extensively used to describe a wide variety of chemical and biological processes [17][18][19][20].
The existing model in [21] has been taken as the base model of this paper.This model is characterized by the use of differential equations that play a decisive role in predicting the tissue regeneration process.The authors propose a 3D model that describes the dynamics of stem cells and nutrients in a regenerative therapy of cardiac tissue through two coupled partial differential equations.They analyze the performance of the model by carrying out simulations that demonstrate how the supply of new nutrients modifies the dynamics of stem cells.
The mathematical model described in [21] , which is still a pure-mathematical model, is primarily focused on the description of the evolution of the regenerated tissue geometry and the reduction of the necrotic zone.Indeed, the mathematical scheme is a free boundary problem, where the free boundary represents the border of the necrotic core: the presence of stem cells promotes the regeneration of heart tissue so that the necrotic zone reduces its size.This fact is not easy to handle even from a theoretical point of view.Regardless, from the point of view of the application, the importance of the model lies in the numerical simulations that allow an approximation to the phenomenon, yielding the expected behavior.
Similar to the model in [21], the domain where the diffusion phenomenon occurs is a region of the heart affected by inflammation, represented as D(t).This domain will be modified as regeneration progresses.The first differential problem is a diffusion-convectionreaction equation that describes the dynamics of stem cells, where s(⃗ x, t) represents the concentration of stem cells in D over the interval [0, T].The second problem of the base model describes the dynamics of nutrients n(⃗ x, t).
In this paper, we add the dynamics of chemokines, describing their chemotactic effects on the stem cells.The study is mainly focused on mesenchymal stem cells and nutrient and chemokine concentrations, though it is well-known that the regeneration process is more complex and involves more biological agents, which could be included in the mathematical modeling in the future.
Stem cells have the capacity to regenerate cardiac tissue and improve the function of the damaged heart; indeed, they can divide and become specialized cells, such as cardiomyocytes.Thus, clearly, the impact of the chemokines in the attraction of these cells toward the damaged area of the myocardium becomes very important.The chemokines induce the migration of the stem cells and also stimulate the movement of certain types of white blood cells, attracting them to areas of inflammation to help the body fight infections, external pathogens and other diseases.Hence, we model the dynamics of chemokines given their role as chemoattractants for stem cells.

Methodology
The model presented in this paper is composed of three partial differential equations.The basic model in [21], consisting of two equations describing the stem cells and nutrient dynamics, is modified to also account for the phenomenon of chemotaxis.
Cytokines and adhesion molecules are biological mediators that lead implanted stem cells from blood vessels to migrate to endothelial cells of a tissue.Once the stem cells reach the tissue, they begin to fulfill their function.Chemokines are particular cytokines produced by certain immune cells and other cells in the body, endowed with chemotactic activity.In the event of a myocardial injury, chemokines are secreted by tissue cells in response to damage caused by a heart attack.
After an acute myocardial infarction, the ventricular remodeling stage begins [22,23].One of its causes is the increase in interstitial collagen deposits, and this alteration leads to increased myocardial fibrosis [24].In the meantime, the damaged tissue seeks to recover and the immune system responds by emitting damage signals so that new cells come to its aid.These signals are the chemokines.That is, after the signal is emitted, a cell receives it and begins its relocation to the damaged area, following a gradient given by chemokines.This biological effect induced by chemokines is called chemotaxis.
In the model presented in [21], which was still a simplified model, the radial attraction was introduced in order to describe the movement of the stem cells toward the necrotic core.Actually, as explained above, this attraction is related to the chemotactic effect generated by the chemokines, which are activated in the core and, in part, move to the inflamed area surrounding the infarcted area.Such an effect is directional, and the direction of the induced motion is the one of the steepest increase in the chemokine concentration, i.e., the direction of the gradient.Thus, in developing our model, the simple radial attraction term toward the center of the necrotic core in the first equation of the basic model in [21], following [10], is replaced with a chemotactic term dependent on the gradient of chemokine concentration in the tissue, which influences stem cell movement by guiding them toward the necrotic zone.
Hence, an equation describing the dynamics of chemokines in the tissue is added to the model.The three equations interact with each other in that the dynamics of stem cells are influenced by both the dynamics of nutrient substances and chemokines.Conversely, the consumption of nutrients and deactivation of chemokines are regulated by the presence of stem cells in the tissue and their dynamics.
In this study, we denote concentrations as follows: stem cells as s, nutrients as n and chemokines as q.

Dynamics of the Stem Cells
The mathematical model for the concentration of stem cells is a diffusion-convectionreaction equation given by: The above equation governs the dynamics of stem cells in the spatial domain D, over a temporal interval [0, T].
The domain D, which represents the inflamed zone, is the region where stem cells are transported from an implantation zone to the necrotic core, which is the portion of tissue injured by the infarction.
Tissue regeneration is promoted by the arrival of stem cells at the border Γ of the necrotic zone.As a result, the necrotic core contracts and its border Γ moves inward.In this way, the dynamics of the stem cells are described by a free boundary model with Γ = Γ(t) and D = D(t).
As an initial condition, a zero amount of stem cells is placed in D(t) at the beginning of the phenomenon.As a boundary condition, the number of cells in a portion A of the fixed boundary ∂E is determined by a function BC(t) decreasing in time.We consider the same function introduced in Figure 3 in [21], which is a smoothed decreasing step function, passing from a maximum value s 0 to 0. In this way the concentration of stem cells on A disappears at a rate that can be predetermined in the model.Furthermore, there is no flow at the boundary, except for the portion A, which is indicated by a Neumann condition through a unit normal vector ν pointing outward of D(t).It is expected that when the stem cells migrate and reach N, substances that regenerate the tissue are activated and engraftment occurs in the damaged tissue.Consequently, the cells disappear upon reaching Γ(t), implying a contraction of the necrotic zone.So the border Γ(t) moves inward with a velocity V(⃗ x, t) proportional to the flux of stem cells.Thus, we have the system for stem cell dynamics [21]: • s: stem cell concentration;

Dynamics of Nutrients
The mathematical model for nutrient dynamics is a diffusion-reaction equation, with variable n(⃗ x, t), which represents the nutrient concentration in D over the interval [0, T].In D(t), the variation in nutrient concentration with respect to time ∂n ∂t is determined by the movement of nutrients represented by a diffusive flow.The reaction terms correspond to the consumption of nutrients, that is, the interaction of nutrients with the stem cells.The consumption ratio is determined by a Hill function of exponent 4 [25].
The initial condition describes an amount n 0 of nutrient in D(t) at the beginning of the phenomenon, that is, at t = 0. To indicate that no nutrient leaves or enters the necrotic core, we have a Neumann condition at the border Γ(t).Furthermore, the flow at the boundary ∂E will be given by a function that can be zero if there is no flow of nutrient (in this case a null Neumann condition), or different from zero if there is exit or entry of nutrient across ∂E, as the case may be.In this model we impose the former condition, supposing that the stem cells feed only on the nutrient already present in the domain.In future research, we can consider a nonzero nutrient flow, which can allow us to model the phenomenon of angiogenesis due to the regeneration of the cardiac tissue.We have:

Dynamics of Chemokines
The diffusion of chemokines occurs from the necrotic area into the inflamed area around the infarction; naturally, inflammation is important for the emission of warning signals that indicate the presence of an injury that the body must repair.Chemokines reach receptors on cells and then induce chemotaxis.The variation in the concentration of chemokines over time is determined by their movement in the inflamed area, and their consumption during chemotaxis.
The mathematical model that describes the dynamics of chemokines is a diffusion equation, with variable q(⃗ x, t), which represents the concentration of chemokines in the inflamed area surrounding the necrotic core.
In D(t), according to standard considerations, we determine the variation in the chemokine concentration over an arbitrary control volume Ω with respect to time, ∂ ∂t Ω q dx by the diffusive movement of chemokines as well as by the interaction of chemokines with stem cells, whose rate is modeled by the mass action principle, and the consumption or death rate of chemokines themselves, proportional to the chemokine concentration: The initial condition for chemokines is given by a concentration q 0 in D(t) at the beginning of the phenomenon, that is, at t = 0. To indicate that chemokines do not leave or enter E, we have a Neumann condition on the boundary ∂E.Furthermore, the flux across the boundary Γ(t) will be given by a constant Φ.We have the following mathematical model: where: • q: chemokine concentration in D × (0, T); • s: concentration of stem cells in D × (0, T); • R sq : chemokine consumption coefficient; • R q : chemokine mortality coefficient; • D q : diffusion coefficient of q; • q 0 : initial concentration; • Φ: constant chemokine flux parameter at the border of the necrotic core.

Coupling
The equation for chemokine dynamics has been coupled to the base model, in the domain D = E\N, where E and N are two concentric cubes described in Figure 1.

Mimetic Differences Method for Partial Differential Equations
In this paper, given that our model includes a new variable, a new partial differential equation and its respective boundary conditions, we opted for high-order mimetic methods, based on finite differences, which are easy to implement.Our proposal was to show that, with mimetic differences, it is possible to find a numerical solution to our coupled system of partial differential equations.
The mimetic difference methods that we used were initially described by Castillo and Grone in the article [26], where the mimetic schemes that the authors construct with second-order difference approximation preserve the properties of the continuous operators, with the same order of accuracy in the interior and on the boundary; this makes it possible for the discretized equations to imitate the physical properties of the problem such as the conservation laws.The authors show how to construct high-order approximations for the derivatives using finite volumes on a uniform one-dimensional mesh.The mimetic scheme uses discretizations that have properties analogous to the divergence theorem, local conservation law and global conservation law.
In [27], the operators given in [26] are improved and high-order divergence and gradient operators in dimensions two and three are constructed, using mimetic finite differences on a staggered grid.These operators have been implemented in MATLAB by the authors Corbino and Castillo [27], making available a library with the operators in matrix form called MOLE (Mimetic Operators Library Enhanced).
To obtain the results shown in this paper, we built the program in MATLAB, which allows the use of mimetic operators such as divergence and gradient.To implement mimetic differences, we perform the discretization of the domain, the discretization of the variables with their initial and boundary conditions and the discretization of the equations, using mimetic differences and mimetic operators.
Let us underline that one of the contributions of this article is the use of a new library of mimetic differences called MOLE.In this library it is possible to find the matrix representations of the Corbino-Castillo high-order operators; Divergence, Gradient and Laplacian.Although these operators can be used in curvilinear domains, the problem in using them lies in the implementation of the boundary conditions for complex domains where there are mobile interfaces, such as the border of a necrotic area.In addition, the current version of MOLE is still being revised for non-Cartesian cases.This is the reason why in the model developed in this paper and in the simulations, the domain has been considered as in Figure 1, with a cube E, which contains an initially cubic necrotic core N and an inflamed zone D = E \ N, while the position and the shape of the implantation zone A of the stem cells were varied.

Mimetic Representation of the Model
Following [28], let S, N and Q be the discrete variables that approximate s, n and q in the centers of the elements of the mesh.We obtain the vectorizations of S, N and Q, where D, G and L are the divergence, gradient and Laplacian operators, respectively, in mimetic differences.The mimetic difference schemes for the system of equations are: where the operator diag corresponds to placing a vector on the diagonal of a matrix and obtaining a diagonal matrix.Using one-step discretization for the time derivative, and linearizing the nonlinear part implicitly, the discrete numerical scheme is as follows: where:

Implementation and Results
Given the discretization introduced within the mimetic framework, a numerical method has been implemented in MATLAB, dwelling upon the MOLE library.
In the presentation, the images show the graph with the upper and lower sides inverted, thus the implantation area is shown below, and the dynamics develop upwards.
The placement of two different values of chemokine concentration in each half of the necrotic core was implemented.
Since the geometry for the boundary available within the library allows for Cartesian grids, the geometry of the necrotic core has been modeled as a cube.However, the geometry and the position of the stem cell implantation zone have been varied, to see how in the different cases, the dynamics of the stem cells, moving toward the necrotic core, would result in a regression of the necrotic core due to the effect of stem cells attracted by chemokines.
Moreover, the implantation of stem cells was carried out also on one side of the boundary, as in Figure 1, to avoid symmetric biases in the execution of the method.In the figures presented in this paper, we chose an off-center square (Figures 2-5) and an off-center disk (Figures 6-8).The side of the square and the diameter of the disk of stem cells are about one-third of the side of the cube where the dynamics evolve: the cube contains at its center the necrotic core, which is initially cubic, too, and around it, the chemokines, the stem cells and the nutrient evolve according to the previous equations.
The initial data and parameters used in the numerical experiments are given in Table 1.
Table 1.Initial data and parameters used in the numerical experiments.

Parameter
Value Unit Let us remark that we do not claim their values to be realistic; instead, they were chosen in order to highlight the effects of stem cell therapy on the necrotic zone.
The chemokines flow from Γ(t), and their interactions with the stem cells take place in D(t).
The dynamics of the stem cells clearly show a qualitative movement in the concentration toward the necrotic core: due to the chemotaxis produced by the chemokines, the stem cells move from the implantation zone A, their interactions with the chemokines take place in D(t), then they reach Γ(t) and, after their engraftment, they reduce the necrotic zone, regenerating the cardiac tissue.Moreover, the moving boundary feature of the numerical method shows how the regression surface of the necrotic core evolves as long as the concentration of stem cells grows outside the moving boundary.The effect is that the necrotic core is slowly eroded from the direction from which the stem cells come, of course-this is evident from the picture of the chemokine concentration compared to the stem cell diffusion, after a while, as in Figures 2-5, where A is an off-center square, and in Figures 6-8, where A is an off-center disk.The moving boundary evolves toward the inner side of the necrotic core, from the side from which stem cells diffuse: in Figure 8, we plot this at the same time as that at which the previous pictures were taken.A is an off-center circle.

Discussion
The model here presented, formed by a system of three coupled diffusion equations, describing the evolution of stem cells, chemokines and nutrients, is qualitatively coherent with the expected experimental outputs, related to stem cell therapy, in particular, for what concerns the regeneration of the damaged area.
In this work, we introduced the dynamics of chemokines, restricting the chemotaxis effects mainly to the concentration of mesenchymal stem cells, nutrient concentration and chemokine concentration, Let us remark that the model proposed in this paper corrects in some sense the previous one, described in [21], clarifying in particular, from a physiological point of view, the attraction of the stem cells toward the necrotic area.Thus, as in [21], it can be considered as the proposal of a new model, corroborated by numerical results, waiting for experimental validation.In fact, as already observed, the initial data and parameter values were chosen in order to highlight the effects of the stem cell therapy on the necrotic zone and do not claim to be realistic.Thus, the model must be still considered a purely mathematical model.In the future, the validation of the model, by means of experimental data, will be necessary, with the goal of making the model more realistic.Once we have validated it, we will have the possibility of improving it, building a more general model, which will take into consideration the paracrine effects of the stem cells, too, having clearly in mind that the regeneration process is very complex and involves more biological agents, which can be included in the future, in the mathematical model.Moreover, for a cell biology model to be as realistic as possible, other important phenomena must be added to the modeling, such as the vascularization and the behavior of nutrients other than oxygen.For what concerns vascularization, in this model, as already remarked, we supposed that there is no flow of nutrient at the boundary ∂E, which implies that the stem cells feed only on the nutrient already present in the domain.In future research, we will consider a nonzero nutrient flow, which can allow us to model the phenomenon of angiogenesis, due to the regeneration of the cardiac tissue.
From a theoretical point of view, the study of very important issues, such as the existence and uniqueness of the solution, the stability analysis and the non-negativity of the model, will be considered in future research.
From a numerical point of view, the main contribution of this work is the implementation of mimetic operators for the solution of the model here proposed through the use of the new library MOLE.This is absolutely novel, thus, at this moment, we have no benchmarks to compare against.
In the future, it is our intention to study the reproduction number of the cardiac cells of our model, whose definition and computation, for PDE models such as reaction-diffusion systems, though intriguing, are typically more complicated.

Figure 1 .
Figure 1.D = E\N (white color) is the inflamed zone, where N (blue color), the necrotic zone, is a concentric cube with faces parallel to the faces of the cube E. A (yellow color) is the implantation area of the stem cells in ∂E.

Figure 2 .
Figure 2. (a) Chemokine concentration at t = 750 s; the highest concentration is on the side closest to the stem cells; (b) chemokine concentration at t = 4104 s; the necrotic core regenerates by 3.5%, and the chemokine concentration has decreased in the regenerated part.A is an off-center square.

Figure 3 .Figure 4 .Figure 5 .Figure 6 .
Figure 3. Profile of the concentration of chemokines after 5000 s.A is an off-center square.The reduction in the necrotic core in the lower boundary is evident.

Figure 7 .
Figure 7. (a) Concentration of stem cells after 1623 s; (b) concentration of chemokines after 1623 s.A is an off-center circle.

Figure 8 .
Figure 8. Deformation of the lower boundary of the necrotic core after 1623 s, due to stem cell action.A is an off-center circle.