A Mathematical Model of Regenerative Axon Growing along Glial Scar after Spinal Cord Injury

A major factor in the failure of central nervous system (CNS) axon regeneration is the formation of glial scar after the injury of CNS. Glial scar generates a dense barrier which the regenerative axons cannot easily pass through or by. In this paper, a mathematical model was established to explore how the regenerative axons grow along the surface of glial scar or bypass the glial scar. This mathematical model was constructed based on the spinal cord injury (SCI) repair experiments by transplanting Schwann cells as bridge over the glial scar. The Lattice Boltzmann Method (LBM) was used in this model for three-dimensional numerical simulation. The advantage of this model is that it provides a parallel and easily implemented algorithm and has the capability of handling complicated boundaries. Using the simulated data, two significant conclusions were made in this study: (1) the levels of inhibitory factors on the surface of the glial scar are the main factors affecting axon elongation and (2) when the inhibitory factor levels on the surface of the glial scar remain constant, the longitudinal size of the glial scar has greater influence on the average rate of axon growth than the transverse size. These results will provide theoretical guidance and reference for researchers to design efficient experiments.


Introduction
Spinal cord injury (SCI) is the damage to the spinal cord that results in a loss of function such as mobility or feeling. An injured spinal cord has a poor intrinsic capacity for regeneration, although some functional recovery does occur. The failure to regenerate is caused by a combination of factors, including neuroinflammation, axonal disruption, death of neurons, glial scar formation, the release of myelin-associated inhibitory molecules, and the lack of growth promoting molecules. It is clear that effective treatment will require a multifaceted combination of strategies. A lot of effort has been made to promote CNS regeneration but with only limited success [1][2][3][4][5][6].
Trauma or disease of a nerve in a mature mammal may result in a massive multiplication of glial cells around the damaged region, which eventually form a dense scar. This glial scar plays a dual role as chemical and mechanical barriers to the axonal regeneration of injured neurons [7][8][9][10]. Many chemical factors from glial scars have been identified, including Nogo-A, Nogo-B, Nogo-C, myelin-associated glycoprotein, and oligodendrocyte-myelin glycoprotein, to inhibit axonal growth [11][12][13]. On the other hand, the size, shape, and hardness of the glial scars represent the mechanical aspects that resist the extension of the regenerative axons.
Some experimental strategies have been employed to improve the CNS regeneration. These strategies include (1) complementing favorable promoting factors in CNS (such as nerve growth factor (NGF)) [14]; (2) decreasing expression of the inhibiting factors such as Nogo-neutralizing antibody (IN-1) and chondroitinase ABC (Chase-ABC) [15,16]; and (3) repairing or reconstructing myelin. For example, transplant of Schwann cell into the SCI site provides a permissive cellular substrate which may enable axons to pass through the scar area and find their targets [17]. However, experimental studies to examine many factors simultaneously are often timeconsuming, costly, and laborious. The optimal concentration ratios of various factors for growth of regenerating axons in a suitable microenvironment are valuable information for researchers in experimental design. The technologies of mathematical modeling could provide an opportunity to elucidate the ratio and distribution law of various impact factors.

Computational and Mathematical Methods in Medicine
Several mathematical models have been developed to describe axonal growth [18][19][20][21][22][23][24][25]. These models generally consist of two parts: (1) the reaction-diffusion equation describes the transmission of nerve factors and other guidance molecules during development and (2) the axonal "growth equation" (based on the cell chemotaxis principle) describes the growth of axons determined by the concentration gradient of the guidance molecules. By allowing for noise in axonal guidance cues and randomized changes in axon growth substrates, a stochastic component has been included in the growth equation [19]. A previous study showed a twodimensional finite difference in the solution and calculation program of "parabolic equations with a gradient term" [23]. Another study obtained a large-scale, two-dimensional simulation result using parallel computing [20]. One of our previous studies used a three-dimensional finite difference method and reported simulation result [22]. The Lattice Boltzmann Method (LBM), a numerical method commonly used for fluid simulation [26][27][28][29], directly describes problems and is convenient for implementing parallel computing. Nevertheless, these studies did not consider regenerating axonal growth in a deprived environment.
In this study a numerical simulation method, which was based on the experiment for the SCI repair by Schwann cell transplantation [17], was used to explore the regenerative axons growing along the glial scar surface. In this model, glial scar was simplified to a rotating ellipsoid in which circular target cells and axons were arranged around the upper and lower part, respectively (see Figure 1). The release rate of the regeneration inhibitors on the scar surface and the neurotrophic factors (NTFs) at the target and the size and the shape of the glial scars were controlled in the simulation, respectively. Concentration gradients of three types of diffusible molecules were tested in this model: (1) Type 1 factors, that is, attractive molecules that are released by the target tissues after nerve injury (such as neurotrophic factor-1 and NGF); (2) Type 2 factors, that is, chondroitin sulfate proteoglycans (CSPGs) that are not neutralized by the Chase-ABC on the surface of glial scar; and (3) Type 3 factors, that is, a variety of growth promoting molecules produced by the transplanted Schwann cells on the surface of glial scar (such as laminin, fibronectin, and neural cell adhesion molecules). Type 1 factors played a leading role in axonal regeneration. Types 2 and 3 factors exhibited balanced and coordinated effects. In addition, cross talk between Type 1 factors and Type 2/3 factors took place through signal transduction [30]. Concentrations of all three types of molecules were denoted as 1 , 2 , and 3 . Coupled reactiondiffusion equations represented concentrations of these three types of molecules varying in space and time, in which they were key parameters coaffecting the growth rate and whether regenerating axons grew or stopped elongating. The "growth equation" of regenerating axons was designed using gradient parameters according to the cell chemotaxis principle. Using a mouse model of spinal cord transaction [4,5,17,31], a boundary condition was established. Numerical simulations were performed using the three-dimensional LBM to determine the quantitative relationship between growth velocity of regenerating axons and concentrations of promoters and inhibitors in a deprived environment. This study provides a theoretical reference for designing the related experiments.

Evolution Equations for the Concentration of Promoting and Inhibiting Factors in Microenvironment.
From the physical and mathematical point of view, three types of diffusible molecules after SCI, the same as nervous system in developmental stages, should be subjected to the first law of Fick [18,[20][21][22]. The equations are Expressions (1)-(3) are multicomponent nonstationary reaction-diffusion equations with nonlinear coupling point sources, where ∇ 2 = 2 / 2 + 2 / 2 + 2 / 2 (the Laplace operator) and r = i+ j+ k ( , , and are coordinates, and i, j, and k are unit vectors in Cartesian coordinate system). 1 , 2 , and 3 , the function for diffusible molecules at position r ( m) at time (s), are the concentration of Types 1, 2, and 3 of diffusible molecules, respectively. 1 , 2 , and 3 are corresponding diffusion coefficients (constants, m 2 s −1 ). −1 , −2 , and −3 are the linear attenuation coefficients (constants, s −1 ). All ∑s are treated as point sources, and (r) is the Dirac delta function. and are the number of the target cells and axons, respectively. r (stabilizing) is the position of the th target cell, and r (growing with the time) is the position of the th growth cone. Constant 1 ( Ms −1 ) is the release rate of Type 1 factors. 2 ( 1 ) and 3 ( 1 ) are the release rate of Type 2 and Type 3 factors, respectively. These two parameters were expressed as a function of 1 (usually nonlinear) in several specific forms to reflect their interactions with Type 1 factors. 20 and 30 are the basic release rates of Type 2 and Type 3 factors, respectively. Consider 2 = 20 (1 − ) and 3 = 30 , = 1 /( + 1 ). , a dimensionless quantity, is the associativity formula of receptor-ligand.
and (1 − ) reflect the competitive relationship of succeeding each other between enhancing and inhibitory factors.
is the dissociation constant.

Equations for Axonal Growth.
The movement of an axonal growth cone is chemotactic and biased towards (attractive chemotactic) or away (repulsive chemotactic) from the chemical source [32,33]. The attractive or repulsive action on a growth cone is referred to as a chemotactic force, which is proportional to the gradient of diffusible molecules from the chemical source. If a growth cone is regarded as a particle, the growth rate of the axon depends on the velocity of the particle. Since axonal growth is very slow (∼0.01-0.05 ms −1 ) [32,33], acceleration or inertia forces can be negligible. Therefore, the velocity of the growth cone is directly proportional to the chemotactic force: where r (= i + j + k) is the position ( m) of the th growth cone of axons at time (s), is the number of regenerative axons, is the viscosity coefficient (Pa⋅s), F (= i + j + k) is the total chemotactic force (pN) from all three types of molecules acting on the th growth cone at time (s), is the number of molecular types, and the variables and parameters with subscript relate to the th type of molecules. We defined a dimensionless vector p with proportionality constant to express the chemotactic force, in which is the concentration of type molecules in position r at time , ∇ is the gradient of (∇ = i / + j / + k / , the Hamilton operator), |Δr | = √Δ 2 + Δ 2 + Δ 2 is the scalar quantity of difference of r , and Σ is the sum of over all three types of molecules. The scalar quantity of p can be simplified into Δ / Σ , the relative difference of , and Δ is the absolute difference of across a distance |Δr |. In practice, one can take |Δr | ∼ 20 m as an average diameter of growth cone. Therefore, a bridge between gradient and relative difference of concentration is established through p mathematically, and the corresponding proportionality constant is given a clear physical meaning whose dimension is equal to [force]/[length]. Using this model, the rate distortion of the growth cone is greatly reduced when the growth cone is close to the target cells. In addition, for one-dimensional single-component problems [34][35][36], p can be reduced to = ( / )⋅(Δ / ) or = Δ / which is consistent with the definition of gradient of concentration in some biophysical areas.

The Curved Surface Equation and Condition of Glial Scar
Surface. Neural cell adhesion molecule L1 can be synthesized by Schwann cells [37], and this is considered to be one of the major reasons why axons can grow firmly attached to the surface of glial scars. In order to reflect this performance of Schwann cells, a constraint equation, the curved surface equation of glial scar surface, was included in the present model.
Due to the hypothesis that glial scar is simplified to a rotating ellipsoid, where the horizontal rotation radius = √ 2 + 2 ( , , and are coordinates in Cartesian system), the lengths of horizontal half axle and vertical half axle are and , respectively, and the sphere center coordinate is ( 0 , 0 , 0 ), and the constraint equation of scar surface can be expressed as Thus, the movement of regenerating axons' growth cones is transformed into a particle motion along the surface of the rotation ellipsoid. The constraint relations for the velocity can be obtained by taking the time derivative of both sides of (5); we can get the constraint relations for the velocity: where is the drawing speed, a longitudinal component of velocity of an axon growing along the scar surface, and can be used to represent a growth rate of axons. Using growth equations (4), in th axon can be rewritten as The meanings of symbols in these equations can be found in Section 2.2.

Development of the Mathematical Model and Numerical
Methods. This paper created a mathematical model which included (1)- (7). Based on the experiment for the SCI repair by transplanting Schwann cells [17], a numerical simulation method was used to explore the regenerative axons growing along the glial scar surface.
Equations (2) and (3) revealed that the source terms nonlinearly related to (4) through point r , and the point source was moving. So (1)-(4) were a set of coupled nonlinear partial differential equations which could only be solved numerically.
The solution of this model contained three steps and three methods: the first step, using the LBM [26][27][28][29] to solve for the concentration field of various factors determined by the reaction-diffusion equations (1)-(3); the second step, using the central difference method to solve for gradients of various factors surrounding the growth cone, and axon growth rates would be solved when the solution for the gradient was then inserted into (4); and, finally, using Euler's method to numerically integrate (4) and solving for the axonal growth path.
The simulation was implemented in two ways: one was to change the release rate of inhibitory factors, and the other was to change the scar diameter. Axonal growth rates were simulated and recorded in simulations. Whether regenerating axons could grow across or around the scar tissue to connect with the target cells was also tested.
In nerve cells, there were reliable data in the order of magnitude obtained from in vitro experiments as follows [34][35][36]: the width of the growth cone 10-20 m, the rate of axon growth 0.01-0.05 ms −1 , the diffusion coefficient of NGF 1 ≈ 100 m 2 s −1 , and dissociation constant ≈ 1 nM in which NGF bind with receptors of the growth cone membrane. In addition, the concentration range was at 0.01K d -10K d and minimum relative concentration difference was 1% which would work on the growth cone. However, there were many data which had not been recognized, including point source release rates of NGF 1 , attenuation coefficient −1 , constant growth rate of axon , and viscosity coefficient . For the above reason, the principle that the diffusion velocity of NGF −1 √ 1 / −1 should be greater than the growth cone velocity was used to calculate attenuation coefficient −1 . Then, using the numerical method, how much point source release rate was needed to generate a concentration range of 0.01K d -10K d and a minimum relative concentration difference of 1% could be calculated. How much / was 0.0045 0 needed to achieve an axon growth of 0.01-0.05 ms −1 in this concentration field could also be extrapolated. Parameters of Type 2 and Type 3 factors were given by the relative ratio of NGF.

Influence of Inhibitor Release Rate on Axonal Regeneration.
In this section, the following hypothesis was tested through numerical simulation: although glial scar exists, it does not affect the growth of regenerating axons as long as the inhibitor level on the surface of the glial scar is below a threshold.
In the simulation, the size and shape of the glial scar were set to a fixed value and inhibitor release rates from their surface were varied. Growth rates of regenerating axons and the concentration of influencing factors near growth cones were recorded to quantitatively analyze the relationships between them. To simplify the process, a cubical compartment with a side length of = 6720 m (see Figure 1 The ratios of the basic release rate of Type 2 and Type 3 factors (release ratio) were 2 = 20 / 1 , 3 = 30 / 1 , and 3 = 3%. The values of these parameter were chosen based on the previous work [34][35][36] and calculations discussed in Section 2.4.
From Figure 2 and Table 1, it is clear that the average drawing speed of axons ( defined by (7)) decreased gradually with the increase of the release rate of growth inhibitory factors. When 2 was 3%, all axons were successfully connected to their target cells (see Figure 1(a)). So when 2 was less than 3%, inhibitory factor did not affect axonal regeneration along the glial scar. When 2 was at 4%, 5%, or 6%, the number of axons successfully connected to the target cell gradually reduced and the average drawing speed of axons slowed significantly (see Figure 1(b)). No axons could regenerate when 2 was equal to 7% (see Figure 1(c) and Table 1). These results showed that Type 2 factor was increased gradually and making a greater impact to hinder axonal regeneration with the increase value of 2 when diffusion coefficients, dissipation coefficients, and other parameters remained unchanged.

Influence of Size and Shape of Glial Scar on Axonal
Regeneration. In this section, axon growth behavior along glial scars with different sizes and shapes was explored. For illustration purpose, five different glial scars (five cases), which were named number 1, number 2, number 3, number 4, and number 5, were used. Here, number 1 is horizontal halfaxle (defined in Section 2. In all cases in this section, the ratios of the basic release rate of Type 2 and Type 3 to Type 1 factors (release ratio) were 2 = 20 / 1 = 3 = 30 / 1 = 3%. Figure 3 shows the longitudinal sections for threedimensional concentration field of Type 1 factors when the axons connected with target cells. It is clear from Figures  3(a), 3(b), and 3(e) that the bypassing glial scar to spread remote areas by Type 1 factors becomes much more difficult with the increase of the vertical size of the glial scar. Figures  3(b), 3(c), and 3(d) show the result of a simulation in which the lateral diffusing capacity for diffusion factors has been weakened with the increase of the transverse size of the glial scar, and the concentration-diffusion is also weakened in the whole computational domain. In comparison to the impact of changes in the vertical radius, changes in the transverse size of glial scar have weaker effects on diffusion of these factors. Figure 4 shows the axonal growth rates varying with size and shape of glial scars. It is clear from Figure 4(a) that the average drawing speed of regenerating axons became slower as the size of glial scar increased. The increase in transverse sizes of the glial scar did not cause significant reduction of axon growth rate. However, the increase in longitudinal size of the glial scar caused significant reduction of axon growth rate. For Figure 4(b), in these regenerating axons, one representative was selected and the time required for the representative axon to reach the target domain was calculated; the resulting curves also show that the longitudinal size of a glial scar had greater influence on the axon growth rate or the time required for an axon to reach its target domain.

Conclusions
This study focused on the physical point of view to explore the reasons for inhibition of CNS axon regeneration from the external microenvironment of the nerve cells. Based on the experiment for the SCI repair by Schwann cells transplantation and certain hypotheses, a new mathematical model was built. This model has two main control parameters: (1) the diameter of a glial scar and (2) the ratio of release rate of axon regeneration inhibitors on the scar surface to release rate of NTF from the target cells 2 . From numerical calculations, simulations, and analyses, the model yielded the following outcomes: (1) Regenerating axons could successfully navigate across the glial scar and connect to their target cells with the support of Schwann cells and the guidance of NTFs concentration gradients of the target cells when 2 was less than 3%. (2) Changes in the longitudinal size of the glial scar had greater influence on the average rate of axon growth and the required time for axon to reach its target domain compared to the transverse size of the glial scar.
Implantation of Schwann cells and chondroitinase is now considered to be one of the most promising treatment strategies for SCI repair [31]. Existing experimental results show that the regenerating axons can grow along the glial scar to reach the target cells occasionally [17]. Simulation results in this study are in satisfactory agreement with existing experimental data in certain conditions. We would like to point out that, as the initial stage in developing this model, this study did not take into account the sprouting mechanisms after neuronal injury, the polymerization of the cytoskeletal protein within the growth cone of the regenerative axon, and other internal factors. Additional experimental data could be integrated to improve this model in the future.
The LBM for a three-dimensional numerical simulation was adopted in this study. This new method has several advantages, especially in dealing with complex boundaries, incorporating of microscopic interactions, and parallelization of the algorithm. First, there is no problem in principle if the impact factors are divided into more components and the glial scars with more complex shapes are taken into account. Although using only mathematical method does not discover and identify what is a promoting factor and inhibitory factor, it provides an important complement to experimental work to elucidate the ratio and distribution law of the impact factors. Second, this model provides a means of integrating data obtained from different experiments and laboratories. As illustrated in the last few sections of Results and Discussion, mathematical models can also have considerable predictive capabilities. After a model has been developed and carefully validated, it can be used to predict the inhibitory factor levels and the size of the glial scar that allow regenerating axons to successfully navigate across the glial scar. Ideally, experimental work and model development should be carried out in close association, using the mathematical model to guide the design and to evaluate experiments and using experimental results to improve the model development.