Glial scar size, inhibitor concentration, and growth of regenerating axons after spinal cord transection.

A mathematical model has been formulated in accordance with cell chemotaxis and relevant experimental data. A three-dimensional lattice Boltzmann method was used for numerical simulation. The present study observed the effects of glial scar size and inhibitor concentration on regenerative axonal growth following spinal cord transection. The simulation test comprised two parts: (1) when release rates of growth inhibitor and promoter were constant, the effects of glial scar size on axonal growth rate were analyzed, and concentrations of inhibitor and promoters located at the moving growth cones were recorded. (2) When the glial scar size was constant, the effects of inhibitor and promoter release rates on axonal growth rate were analyzed, and inhibitor and promoter concentrations at the moving growth cones were recorded. Results demonstrated that (1) a larger glial scar and a higher release rate of inhibitor resulted in a reduced axonal growth rate. (2) The axonal growth rate depended on the ratio of inhibitor to promoter concentrations at the growth cones. When the average ratio was < 1.5, regenerating axons were able to grow and successfully contact target cells.

Research Highlights

Larger glial scars and higher inhibitor release rate result in slower axonal growth velocity.

Axonal growth velocity depends on the ratio between inhibitor and promoter concentrations at the position of growth cones. With a ratio < 1.5, regenerating axons were able to grow and successfully connected to target cells.

INTRODUCTION

Spinal cord transection has been shown to induce active gliosis at the wound site and cause a glial scar, which acts as an obstacle and affects neuronal connections at both ends[123456]. Many factors have been shown to inhibit axonal regeneration[789], such as Nogo-A, Nogo-B, Nogo-C, myelin-associated glycoprotein, and oligodendrocyte myelin glycoprotein. Numerous studies have demonstrated that levels of chondroitin sulfate proteoglycan family members are upregulated following central nervous system injury[51011]. In addition, some drugs and techniques that have been employed to block these inhibitors and promote growth of regenerating axons, such as Nogo-neutralizing antibody (IN-1)[1213], chondroitinase ABC[11], and genetic deletion[14]. If promoters (such as laminin) are simultaneously removed, growth of regenerating axons is also suppressed[515]. Moreover, the addition of abundant neurotrophic factors in a deprived environment can enhance growth of regenerating axons, as shown in a study utilizing transplanted peripheral nerve tissue as the target tissue[516]. However, the promoting or inhibitory effects of a factor on axonal growth is not a simple “all or none” phenomenon[1718]. Each factor functions within a range, and the factors can interfere with each other or “cross-talk” when combined[51519]. Therefore, further studies are needed to determine the concentration ratios of various factors for growth of regenerating axons in a suitable microenvironment. Mathematical techniques could be utilized to determine these ratios.

Several mathematical models have described axonal growth, and these models are primarily composed of two parts[202122232425]. (1) The reaction-diffusion equation describes transmission of nerve factors and other guidance molecules during development. (2) The axonal “growth equation” (based on the cell chemotaxis principle) was determined according to a concentration gradient of these guidance molecules. By allowing for noise in axonal guidance cues and randomized changes in axon growth substrates, a stochastic component has been included to the growth equation[21]. A previous study showed a two-dimensional finite difference in the solution and calculation program of “parabolic equations with a gradient term”[25], and another study obtained a large-scale, two-dimensional simulation result using parallel computing[22]. Our previous study reported a simulation result using the three-dimensional finite difference method[24]. Nevertheless, these studies did not consider regenerating axonal growth in a deprived environment. The Lattice Boltzmann method, a numerical simulation method, has been used to simulate flow[26272829], and the method directly describes problems and is convenient for implementing parallel computing.

The main strategy of this study is as follows: factors that affect axonal regeneration in a deprived environment are divided into three types: (1) type 1, promoters associated with target tissues (such as neurotrophic factor-1 and nerve growth factor); (2) type 2, inhibitor in scar tissue (such as chondroitin sulfate proteoglycans and other extracellular matrix molecules semaphorin 3, ephrin-B2, and Slit proteins); (3) type 3, helper factor in the substrate (such as laminin, fibronectin, and neural cell adhesion molecules). Type 1 factors play a leading role in axonal regeneration. Type 2 and 3 factors exhibit balanced and coordinated effects. In addition, cross-talk between type 1 factors and type 2/3 factors takes place through signal transduction[19]. Concentrations of all three types of molecules have been denoted as ρ1, ρ2, and ρ3. Their temporal and spatial evolution has been described using a group of coupled reaction-diffusion equations. The “growth equation” of regenerating axons was designed using gradient parameters according to the cell chemotaxis principle. Using a mouse model of spinal cord transection[453031], a boundary condition was established. Numerical simulation was performed using the three-dimensional lattice Boltzmann method to determine the quantitative relationship between growth velocity of regenerating axons and concentrations of promoters and inhibitors in a deprived environment.

RESULTS

Uncertain parameter values in the mathematical model

The mathematical model established in the present study and compiled computer programs contained numerous uncertain parameters, which should be identified by “presupposition-calculation-correction.” Some precise data were obtained by in vitro tests[323334]: 10–20 µm growth cone width, 0.01–0.5 μm/s axonal growth speed, diffusion coefficient of nerve growth factor (belongs to type 1 factors identified in the study) D1≈100 μm2/s, dissociation constant of nerve growth factor combined with the growth cone membrane receptor K ≈ 1 nM. In addition, the growth cone-affecting concentration ranged from 0.01–10 K, and the minimal relative concentration difference was 1%. However, the nerve growth factor point-source release rate σ1, absorption coefficient k-1, velocity coefficient λ, and viscosity coefficient µ did not result in accurate data. The value k-1 was calculated according to the principle that nerve growth factor diffusion velocity should be greater than growth cone velocity. It was possible to speculate the point source release rate required in a concentration field where the concentration ranged from 0.01 K–10 K and the minimal relative concentration difference was 1%. It was also possible to speculate on λi/µ results when axonal growth speed reached 0.01–0.5 μm/s in this type of concentration field. The type 2 and 3 factor parameters were calculated according to the nerve growth factor relative ratio. The order of magnitude of the ratio of diffusion coefficient to absorption coefficient was . The order of magnitude of the ratio of velocity coefficient of the axonal growth cone was λ1/λ2 = −λ1/λ3 = −1 λi/µ=1 µm/s. The order of magnitude of the ratio of internal diameter of vertebral canal to the length of obtained specimen was α = d/l = 0.275. The ratio of the diameter of glial scar to internal diameter of vertebral canal (diameter ratio) was β = d/l. The ratio of the release rate of type 2 and 3 factors to the release rate of type 1 factor (release ratio) was η2 = σ20/σ1 and η3 = σ30/σ1, respectively. The original value of the release ratio η2, η3 of type 2 and 3 factors in internal and external tissues of the glial scar was 2%. However, η2, η3 in external tissues and η3 in internal tissues were constant (2%), but η2 in internal tissues increased, resulting in elevated inhibitor concentrations in the scar. In other words, there were only two control parameters: β and η2, which were used to regulate glial scar size and inhibitor concentration. However, the remaining parameters were constant. Based on these results, the numerical calculation of the research content was conducted and results were obtained.

Concentration distribution of positive growth factor and inhibitor, growth course, and growth speed of regenerating axons after spinal cord transection

Figure 1 shows the distribution profiles of the promoter concentrations (ρ1) generated by target cells and the inhibitor concentrations (ρ2) generated by the glial scar. The concentration was high in the bright region; the concentration was low in the dark region; and the concentration outside the vertebral canal was zero. ρ1 and ρ2 exhibited a gradient distribution in the vertebral canal. Gradients of ρ1, ρ2, and ρ3(not shown in Figure 1) determined the course and speed of axonal growth.

Figure 1

Slice of three dimensional fields of the promoter concentration (ρ1) and inhibitor concentration (ρ2) generated by the target and the scar, respectively.

Concentrations are high in the bright region; concentrations are low in the dark region; and concentrations outside the vertebral canal (dark blue region) are zero.

ρ1 and ρ2 in the vertebral canal exhibit a gradient distribution. The coordinate axis X, Y, and Z represent the size of the geometric model of spinal cord transection.

As shown in Figure 2, the sphere represents the glial scar following spinal cord transection; the semi- cylindrical shell represents the vertebral canal. To show growth of regenerating axons, the top portion was removed, but an integral cylindrical shell was utilized during calculations. New axons, which sprouted from remaining neurons, were observed in the left portion of the glial scar, and remaining target cells were found in the right side. Figure 2A shows a blank control. Glial scars existed, but the inhibitor release rates were identical in internal and external tissues (η2 = 2%). At this time, axonal growth was similar to growth observed during nervous system development, with the presence of axonal fasciculation and defasciculation[202124]. Ultimately, the axons reached their target cells. Figures 2B and C show that a small glial scar was not able to block growth of regenerating axons when inhibitor release rates in the glial scar were identical (η2 = 5%). Glial scars with an equal diameter (β = 0.582) did not inhibit growth of regenerating axons when the release rates of inhibitors were small (Figures 2A, B, D).

Figure 2

Influence of glial scar size (β) and enclosed inhibitor release rate (η2) on growth course of regenerating axons.

The sphere represents the glial scar following spinal cord transection; the semi-cylindrical shell represents the vertebral canal. New axons sprouting from remaining neurons are seen in the left side of the glial scar, and remaining target cells are detected in the right side.

No.1 refers to the number of axons. Color shades do not have physical or chemical meanings. Coordinate axis X, Y, and Z represent size (length, width and height) of geometric model of the spinal cord transection.

A: β = 0.582, η2 = 2%; B: β = 0.582, η2 = 5%;

C: β = 0.418, η2 = 5%; D: β = 0.582, η2 = 4%.

A is a blank control. Glial scars exist, but the release rates of inhibitors are identical in internal and external tissues (η2 = 2%); B and C show that a small glial scar does not easily block growth of regenerating axons when release rates of inhibitors in the glial scar remain unchanged (η2 = 5%).

A, B, and D show that glial scars of an equal diameter (β = 0.582) do not easily block growth of regenerating axons when release rates of inhibitors are small.

Figures 3A, 4A, and 5A list concentrations of type 1–3 factors at the position where axon No. 1 existed in Figure 2A. In addition, the figures show changes in movement speed of the growth cone and time of axonal growth. A in all figures refers to normal axonal growth (blank control), glial scars did not chemically or physically differ from the common substrate). Promoter concentrations (Figure 3A) increased exponentially with time as the growth cone reached the target cells (release source), and concentrations greatly fluctuated near the target cells. Axons reached the target cells, which generated promoters and transported them into the neuronal cell body via the axon. Therefore, promoter concentrations were low around the target cells and, therefore, not attractive for other axons. In addition, other axons grew towards the region with high concentrations of the promoter. Inhibitor concentration (ρ2) slightly differed from the helper factor concentration (ρ3) generated by the substrate (including glial scar), which was determined by the release rate pattern of type 2 and 3 factors (equations 1–3 in the methods section). A balance point appeared at approximately 3 000 minutes (Figure 4A), and an additional balance point appeared when the axons successfully contacted the target cells. Axonal growth velocity (Figure 5A) was determintly determined according to promoter concentration gradient (ρ1), which was generated by the target cells, and velocity changes were consistent with ρ1 changes. At the beginning, target signals were weak, but then slowly advanced and sped up. Influenced by the connection of other axons to target cells, the touch time was increased and forward velocity became slow. Subsequently, velocity increased when the axons nearly reached the target cells, although velocity remained within 0.01–0.5 μm/s. All axons in Figure 2A successfully contacted target cells; this took 4068 minutes.

Figure 3

Changes in promoter concentrations (ρ1, Y-axis), as recorded by tracking the No. 1 growth cone of regenerating axons with the time (X-axis) of axon growth. A and B, respectively, correspond to growth cone of No.1 axon, as shown in Figures 2A and B.

A (Figure 2A) refers to normal axonal growth (blank control). ρ1 increases exponentially in time to when the growth cone reaches the target cells, but fluctuates during the connection.

B (Figure 2B) represents large glial scar and high release rate of inhibitors.

ρ1 gradually increases with time to when the growth cone reaches the glial scars, and then ρ1 maintains no changes and the growth cone ceases to grow.

Figure 4

Changes in inhibitor concentration (ρ2, Y-axis) and helper factor concentration (ρ3, Y-axis) in the No.1 growth cone of regenerating axons with time of axonal growth (X-axis). A and B refer to growth cone of No.1 axon in Figures 2A and B, respectively.

A (Figure 2A) refers to normal axonal growth (blank control). The inhibitor concentration ρ2 and helper factor concentration ρ3 are almost mirror symmetry with time to successful connection of growth cone and target cells. A balance point appears immediately prior to connection and upon successful connection.

B (Figure 2B) represents a large glial scar and high release rate of inhibitors. ρ2 represents a process of increase, fluctuation, increase, decrease, and a stable increase over time to when the growth cone ceases to grow, suggesting that axons exhibit motility and avoid high inhibiting concentrations.

However, the axons stop growing due to high inhibitor concentrations in the microenvironment, and helper factor concentration ρ3 was low.

Figure 5

Changes in growth rate (Y-axis) of the axonal growth cone (No. 1) with time of axonal growth (X-axis).

A and B refer to growth cone of No.1 axon in A and B, respectively. Vc: total velocity; Vcx: forward velocity; Vcy and Vcz: lateral velocity.

A (Figure 2A) refers to normal axonal growth (blank control). At the beginning, the movement is very slow because target signals are weak, but then slowly advances and speeds up.

Influenced by the connection between other axons and target cells, the touch time increases, forward velocity slows, and then velocity increases when the axons are close to the target cells, although velocity remains within 0.01–0.5 μm/s.

B (Figure 2B) shows a large glial scar and high release rate of inhibitors. Forward velocity gradually accelerates over time to when the growth cone reaches the glial scar. The dashed line in B represents lateral movement of the growth cones.

Lateral velocity is significant once the growth cone reaches the glial scars. Finally, the growth cone stops growing, because the growth cone was not able to overcome the scar impermeability.

Figures 3B, 4B, and 5B show concentrations of type 1–3 factors in the position where axon No. 1 existed in (Figure 2B), and velocity of the growth cone varying with time of axonal growth. Axonal growth was immediately recorded after the experiment. Axons reached the glial scar at 1 000 minutes, but stopped growing at 1 500 minutes. The abrupt increase in inhibitor concentration (Figure 4B) and slow growth speed (Figure 5B) represented the appearance of a glial scar. Unchanged promoter (Figure 3B) and inhibitor (Figure 4B) concentrations, as well as zero speed of axonal growth (Figure 5B), were used to represent growth inhibition. Figure 4B shows that inhibitor concentrations at the growth cone underwent increases, decreases, and stable increases once the axons (No. 1) reached the glial scars. These results suggested that axons avoided high concentration of inhibitors and ceased to grow as a result of high inhibitor concentrations in the microenvironment. Over time, the growth (Figure 5B) slowed and then completely ceased.

Relationship of growth velocity of regenerating axons to glial scar size and promoter and inhibitor concentrations surrounding the growth cone

Figure 6 shows the relationship curve of mean growth velocity to the inhibitor release ratio (η2) and glial scar diameter ratio (β) in an axon tract of the image from Figure 2.

Figure 6

The relationship curve of mean growth velocity (Y-axis) to inhibitor release ratio (η2, X-axis) and glial scar diameter ratio (β) in an axon tract of Figure 2.

Axonal growth velocity decreases with increased inhibitor release ratio η2. If the inhibitor release ratio η2 is identical, axon growth velocity is low when the glial scar diameter ratio (β) is large.

Axon growth velocity decreased with increased inhibitor release ratio. If the inhibitor release ratio was unchanged, the inhibitor concentration ratio (ρ2/ρ1)mean increased with the size of the glial scar (Figure 7). Figure 8 shows a scatter diagram of the relationship between mean inhibitor concentration ratio and mean growth velocity of regenerating axons.

Figure 7

The relationship curve of mean inhibitor concentration ratio ((ρ2/ρ1)mean, Y-axis) to inhibitor release ratio (η2, X-axis), and the glial scar diameter ratio (β) in an axon tract of Figure 2.

Inhibitor concentration increases with an increasing release ratio η2. If the inhibitor release ratio η2 remains unchanged, the inhibitor concentration increases with an increasing glial scar diameter ratio (β).

Figure 8

Scatter diagram of the relationship between mean inhibitor concentration ratio ((ρ2/ρ1) mean, X-axis) and mean growth velocity (Y-axis) of regenerating axons in an axon tract of Figure 2.

Regardless of the size (β) of the glial scar or the inhibitor release rate (η2), if (ρ2/ρ1) mean is small, regenerating axons will rapidly grow. In contrast, with an increasing (ρ2/ρ1) mean, regenerating axons slowly grow.

As shown in Figure 8, regardless of the size (β) of glial scar or inhibitor release rate (η2), if (ρ2/ρ1)mean was small, regenerating axons underwent rapid growth. In contrast, if (ρ2/ρ1)mean was large, the regenerating axons would slowly grow. Axons ceased to grow, which was similar to Figures 2B and 5B, where Vmean < 0.005 μm/s and (ρ2/ρ1)mean > 1.5 (shown in lower right corner of Figure 8).

DISCUSSION

The present study did not analyze secondary injury or apoptosis following neuronal injury[353637], sprouting mechanisms of remaining neurons, or polymerization of skeleton protein in growth cone of regenerating axons. Nevertheless, this study investigated the influence of the inhibitory effect of an external microenvironment on axonal regeneration, including concentration and time.

The major difference between this study and previous studies was as follows: variable descriptions in mathematics should be redefined in a deprived environment. The inclusion of a boundary condition (in particular, the existence of glial scars) was complicated[20212223242526272829]. The Lattice Boltzmann method was first utilized in such a mathematical model, and the condition or assumption for the present study was: (1) remaining neurons could survive and sprout new axons by utilizing suitable measures following central nervous system injury. (2) Mechanisms underlying axonal growth and guidance were identical to those during neural development. Axonal growth cones perceive a gradient of guidance molecules and exhibit a pathfinding capacity. (3) Following central nervous system injury, surviving neurons do not structurally contact target tissues, although the target cells secrete neurotrophic factors, which diffuse outwards. (4) Following central nervous system injury, glial scars and inhibitors, which suppress axonal growth in the external microenvironment, can diminish and be degraded by medical interventions.

Results from a previous study confirmed that target tissues provide nutrition to neurons via axons by utilizing reverse transport[35]. However, this supply system becomes damaged following axonal injury. Even if target tissues synthesize and secrete neurotrophic factors, transport efficiency remains low due to diffusion. If specific measures are not taken, the remaining neurons would “starve to death” or become apoptotic. Some measures not only prevent secondary injury, but also contribute to sprouting. For example, olfactory ensheathing cells, which exhibit a lifetime regeneration capacity and are rich in neurotrophic factors, were implanted into injured sites[238]. The second assumption is the physical mechanism of model establishment, which is an indisputable fact during developmental stages[3940]. However, the third assumption likely does not hold. Exogenous neurotrophic factors, such as nerve growth factor and brain-derived neurotrophic factor, have been injected into injured target tissues[516]. Exogenous histiocytes that secrete these neurotrophic factors can be implanted[516]. These histiocytes are guidance molecules and induce axon growing towards target cells during development. The fourth assumption has been supported by previous successful studies, such as myelin-associated neurite growth inhibitory protein antibodies (IN-1)[1213], chondroitinase ABC[511] and genetic deletion[14].

In the present study, molecules affecting axonal regeneration were assigned to three types: (1) promoters associated with target tissues; (2) inhibitors in scar tissues; and (3) helper factors associated with the substrate. However, the actual conditions were complicated, and some molecules represented secretory, diffusible, and transmembrane types[3940]. The effects on axonal regeneration were promoting or inhibitory, with determinacy and randomness.

The obstruction effects of glial scars on regenerating axons are mechanical or chemical, and the present study mainly analyzed the concentration of chemical factors. Mechanical obstruction is embodied by size and firmness of glial scars, which, in the present study, implied an axonal growth velocity coefficient λ and substrate viscosity coefficient μ from a mathematical model. The scars became firm at 9 months following central nervous system injury[1841], and regenerating axons did not grow into the scars. However, previous studies demonstrated that olfactory ensheathing cells permeate into the scars and serve as a substrate. These results suggested that regenerating axons traverse this substrate.

In the present study, the geometries of the vertebral canal and glial scar were simplified to a cylindrical shell and sphere, respectively. Figure 8 illustrates that (ρ2/ρ1)mean affects axonal regeneration; β and ρ2 are a means to regulate (ρ2/ρ1)mean. The shape of the glial scar did not affect regeneration. If (ρ2/ρ1)mean was controlled, axonal regeneration would also be controlled.

The present study used a mouse model of spinal cord transaction. Mathematical model was established using known experimental data and certain assumption. There are two control parameters: (1) the ratio of glial scar diameter to internal diameter of vertebral canal (β); (2) the ratio of inhibitor release rate to neurotrophic factor release rate (η2). There are two evaluation indices: (1) mean growth velocity of regenerating axons Vmean; (2) mean ratio of inhibitor concentration to neurotrophic factor concentration (ρ2/ρ1)mean. Numerical calculation and analysis revealed that: (1) a smaller glial scar size and decreased inhibitor release rate resulted in decreased inhibitor concentration and smooth growth of axons. (2) With a ratio of inhibitor concentration to neurotrophic factor concentration (ρ2/ρ1)mean < 1.5, the axons smoothly grew and reached the target cells.

MATERIALS AND METHODS

Design

Mathematical model study.

Time and setting

Experiments were performed at the Laboratory of Applied Mathematics and Mechanics, Shanghai University, China from December 2010 to June 2011.

Materials

Evolution equation of promoter and inhibitor concentrations following spinal cord transection

The molecules affecting axonal regeneration were characterized by three types: (1) promoters associated with target tissues; (2) inhibitors in scar tissues; and (3) helper factors associated with the substrates. The diffusion should obey Fick's law, and its form is identical to the stage of nervous system development[202122232425]:

Equations (1)–(3) are unstable nonlinear reaction-diffusion equations with moving point sources. ▽2 = ∂2/∂x2 +∂2/∂y2 +∂2/∂z2 is the Laplace operator. ρ1, ρ2, and ρ3 are concentrations (µM) of Type 1, 2, and 3 factors, respectively, and they are functions of spatial position r (µm) and time t(s), and r = xi + yj + zk (i, j, k represent unit vector). D1, D2, and D3 are diffusion coefficients (μm2/s) and are constants. k−1, k−2, and k−3 represent the linear absorption coefficients (s-1) and are constants. The point source is represented by the term Σ in each formula, where δ is the Dirac function, i.e., δ(0) = 1, δ(else) = 0. N is the number of target cells. N is the number of axons. is the position (stationary) of the jth target cells. is the position (changes with time) of the kth growth cone. σ1 is the release rate of Type 1 factor (μM/s) and is a constant. σ2(ρ1) and σ3(ρ1) are release rates of Type 2 and 3 factors, respectively, and are the nonlinear functions of ρ1 : σ2 = σ20(1−R), σ3 = σ30R, and R = ρ1/(K + ρ1), where σ20 and σ30 are the basic release rates of Type 2 and 3 factors (μM/s). R obeys the ligand-receptor binding law. R and (1–R) reflect the increase and decrease relationship between the promoter and inhibitor. Kd is the dissociation constant.

Axonal growth equation

Growth cone movement is a chemotactic process[3940]. Attraction and repulsion[3940] during growth cone movement are represented by force F(r) and are directly proportional to the concentration gradients of the attractants and repellents in microenvironment. If the growth cone was simplified to a particle, growth velocity would depend on stable movement velocity of growth cone. Axonal growth is very slow (0.01 μm/s), and acceleration or inertia could be neglected. Therefore, the velocity of growth cone movement is directly proportional to force:

where µ is the viscosity coefficient (Pa·s) and i is the number of components. There are three types of factors and three components. represents the spatial position of the axonal growth cone k, ||Δr||= , refers to a relative concentration difference, and ||~Δρ/ΔΣ represents the ratio of concentration difference Δρ of the ith types of factors to the sum ΔΣ of concentrations of Type 1–3 factors at the position of the axon growth cone k. ▽ = ∂/∂x + ∂/∂y +∂/∂z represents the Hamilton operator. λ is a proportionality constant (velocity coefficient), and is positive for attraction and negative for repulsion, which specifically needs to be validated by experiments. Numerous studies[202122232425] have directly utilized the mathematical gradient ▽ρ to describe chemotaxis, i.e., , therefore, velocity distortion was observed when the growth cones reached the target cells[202122232425]. The present study utilized to represent chemotaxis, which is sensitive to relative concentration differences (rather than absolute concentration values). Moreover, in the present study, the physical concept of the proportionality constant is clear, and the dimension is [force]/[length], and the velocity distortion has been greatly improved. In the one-component one-dimensional problem[323334], the relationship between relative concentration difference and gradient is p = ∂ρ/ρ=(∂ρ/∂x)·(Δx/ρ).

The present study did not analyze sprouting mechanisms following neuronal injury or polymerization of skeleton protein in growth cones. Future studies are needed to determine the intrinsic factors of axonal regeneration.

Methods

Numerical simulation of the three-dimensional Lattice Boltzmann method

The boundary condition was established in accordance with a mouse model of spinal cord transection[453031]. Geometries of vertebral canal and glial scar were simplified to a cylindrical shell and sphere, respectively. Sphere size is a mechanical property of the scars. As displayed in equations (2) and (3), the source terms are functions of r(x, y, z) and are associated with equation (4). In addition, the point sources are movable. Therefore, equations (1)–(4) represent a system of partial differential equations (coupled, nonlinear), which should have to be solved using a numerical method. The calculation program was compiled using Matlab 7.4 software. Figure 1 was generated utilizing Tecplot 10, in which the data was from the Matlab calculations. The remaining figures were generated following calculations utilizing Matlab 7.4. The computational procedures included three steps, and three numerical methods were used. (1): concentration fields of each factor described by reaction-diffusion equations (1)–(3) were determined using the lattice Boltzmann method[26272829]. (2): the gradient of each factor surrounding the growth cones was determined using the central difference method. Results were substituted in the axonal growth equation (4), resulting in growth velocity for each axon. (3): numerical integration was introduced to solve the equation (4) using the Euler method, resulting in an axonal growth course.

The present study hypothesized that the globular glial scar was located in the vertebral canal center, and regenerating axons were at one side of the scar. In addition, the target cells that secreted neurotrophic factors were at another side of the scar. Inhibitor levels were greater in scar tissues compared with outside the scar. The present study investigated the inhibitor release rate in the scars, as well as the scar diameter.

In the numerical simulation, the geometrical size and time measurement of a model were not necessarily identical to the organism prototype and were scaled to fit the computer. Therefore, equations (1)–(4) were scaled non-dimensionally using suitable length scale L0, time scale T0 and concentration scale M0. Following a dimensionless paradigm, the equations still utilized the original symbols.