Numerical characterization of regenerative axons growing along a spherical multifunctional scaffold after spinal cord injury

Spinal cord injury (SCI) followed by extensive cell loss, inflammation, and scarring, often permanently damages neurological function. Biomaterial scaffolds are promising but currently have limited applicability in SCI because after entering the scaffold, regenerating axons tend to become trapped and rarelyre-enter the host tissue, the reasons for which remain to be completely explored. Here, we propose a mathematical model and computer simulation for characterizing regenerative axons growing along a scaffold following SCI, and how their growth may be guided. The model assumed a solid, spherical, multifunctional, biomaterial scaffold, that would bridge the rostral and caudal stumps of a completely transected spinal cord in a rat model and would guide the rostral regenerative axons toward the caudal tissue. Other assumptions include the whole scaffold being coated with extracellular matrix components, and the caudal area being additionally seeded with chemoattractants. The chemical factors on and around the scaffold were formulated to several coupled variables, and the parameter values were derived fromexisting experimental data. Special attention was given to the effects of coating strength, seeding location, and seeding density, as well as the ramp slope of the scaffold, on axonal regeneration. In numerical simulations, a slimmer scaffold provided a small slope at the entry “on-ramp” area that improved the success rate of axonal regeneration. If success rates are high, an increased number of regenerative axons traverse through the narrow channels, causing congestion and lowering the growth rate. An increase in the number of severed axons (300–12000) did not significantly affect the growth rate, but it reduced the success rate of axonal regeneration. However, an increase in the seeding densities of the complexes on the whole scaffold, and that in the seeding densities of the chemoattractants on the caudal area, improved both the success and growth rates. However, an increase in the density of thecomplexes on the whole scaffold risks an over-eutrophic surface that harms axonal regeneration.Although theoretical predictions are yet to be validated directly by experiments, this theoretical tool can advance the treatment of SCI, and is also applicable to scaffolds with other architectures.


Introduction
Spinal cord injury (SCI) typically results in a permanent loss of neurological function below the level of the injury [1,2]. Recovery from SCI is poor owing to weak intrinsic growth capacity of axons and an inhibitory microenvironment, as well as a lack of suitable growth substrates and growth-stimulating factors, that limit axonal regeneration. Many efforts have been made to overcome these impediments [3,4,5,6,7], including implanting of cells [8,9] and/or biomaterials [10,11,12,13,14,15,16,17], delivery of growth factors and degradation of inhibitory matrix molecules [18,19,20,21,22,23,24,25], activation of an intrinsic growth program [26], and stabilization of growth cones and the axonal cytoskeleton [27]. Biomaterial scaffolds that can fill or bridge a lesion cavity, provide a substrate for cell seeding, offer physical guidance to regenerating axons or act as a vehicle for drug delivery are highly promising for cellular and molecular regenerative therapies [19]. However, the use of biomaterial scaffolds as a bridge for SCI repair is complicated; for instance, a popular biomaterial scaffold has tunnels/linear pores, which could guide regenerating axons along these tunnels [14,19]. However, significant congestion was observed at the entry points to the scaffold, which lowered the number of axons entering the pores [19]. If pores were coated or seeded with extracellular components (ECM) components [14] or cells [19], the axons were likely to enter the tunnels from both rostral and caudal ends, but became trappedwithin them [14], unless an additional injection of cells and/or growth factors close to either side of the bridge was administered [18,19], after which some of the trapped axons would be attracted to the injection site and re-enter host tissues [19]. These issues involve the architecture of the "on-ramp" and "off-ramp" parts of the bridge [3], and the concentration gradients of the molecules released from the cells around and seeded on the scaffold, and their influence on successful growth rates of regenerative axons have not yet been clarified.As an analytical tool, mathematics can be used to address many aspects of these issues, for example, the role of slopes at the on-ramp part of a scaffold. When a scaffold with linear pores is used, the axons whose growth cones face to the pores can easily enter the pores owing to the zero entry slope for them, whereas for axons beside the pores, an abrupt or infinite (mathematically) slope exists, which obstructs the axons and leads to congestion. Therefore, theoretical models which address the physical basis underlying the regulatory effect of ligand gradients on axonal growth cone motility are highly desirable, astheyhave been previously used for studying axonal growth during neural development [28,29,30,31]. Axonal growth or regrowth follow the same principle, both during development and following injury. A difference between them, however, has been observed in the level of factors present in their microenvironments [1,2]. Therefore, a theoretical model for axonal growth during development should be modified for studying SCI.Previously, in an initial study [32], we numerically demonstrated that SCI leaves a spherical glial scar. During therapeutic treatment with Schwann cell coating, regenerating axons were shown to grow across the scar and reach their target cellsunder certain conditions. In this study, we assumed a solid, spherical, multifunctional, biomaterial scaffold that would bridge the rostral and caudal stumps of a completely transected spinal cord in a rat model and guide the regenerative rostral axons toward the caudal tissue (Fig 1). The rostral and caudal areas of the scaffold were referred to as the entry "on-ramp" area and the exit "offramp" area, respectively. The whole scaffold was coated with extracellular matrix components, and the caudal area (off-ramp) was additionally seeded with chemoattractants. In this context, the chemical factors on and around the scaffold were formulated to several coupled variables, and the parameter values were derived fromexisting experimental data. The effects of axonal regeneration on the on-ramp slope of the scaffold and biomedical modifications were numerically simulatedto provide a quantitative interconnection between axonal regeneration and the biophysical properties of the components for rational design of scaffolds in tissue engineering [33], as well as contribute to a better understanding of the biological processes involved. From this model and the simulations, we can learn what constitutes a good scaffold, good entry points for regenerative axons,and the resultant concentrations of the molecules from the microenvironment and the scaffold, and how the sources of the main attractants vary along the length of the scaffold in a monotonic curve. Moreover, our model can be modified to optimize scaffolds with other architectures, includingdifferences in size and shape, in concentration of molecules inside and around the scaffold, and before and after experiments.

Materials and methods
Materials and methods have been described primarily from the computational point of view.

Geometry of the scaffold
The scaffold is assumed to be a rotational ellipsoid in such a coordinate system that the z-axis is the rotational axis and is orthogonal to the x-and y-axes (Fig 1B), and the origin O (0, 0, 0) is at the ellipsoid center. The lateral and longitudinal semi-axes of the ellipsoid are denoted by r a and r b , respectively. When r a and r b are equal, the ellipsoid becomes a perfect sphere. Moreover, the rostral and caudal sides of the scaffold are defined as the entry on-ramp and exit offramp, respectively [3]. The longitudinal length of the on-ramp/off-ramp (measured from the scaffold end-point to where the rotational radius of the scaffold surface equals r R /r C ) is r b /3. The length parameters (r a ,r b , r C ,r R , and others) are scaled to dimensionless values by the characteristic length L of the model. The parameter r a is controllable, and r b is fixed at 2 × r b × L = 3 mm (spanning the gap of the spinal cord). Three scaffolds were prepared: slim (r a ×r b = 0. For brevity, all scaffolds are referred to as "spherical", noting that when the longitudinal size r b is fixed, the on-ramp/off-ramp slope is proportional to the lateral size r a . Correspondingly, the entry slopes of the slim, round, and stocky scaffolds are small, medium, and large, respectively. Therefore, when studying the effect of the entry slope of the scaffold on axonal regeneration, we need to vary r a in the model. The scaffold architecture, coordinate system, and symbols employed in the mathematical models. The fabricated scaffold is assumed to be coated with hydroxylapatite (HA)/extracellular matrix (ECM) components (collagen I, fibronectin, and laminin I) and HA/ LV-chondroitinase ABC (ChABC) over the whole surface (green). Additional HA/LV-NT-3/LV-brain-derived neurotrophic factor (BDNF) was coated on the off-ramp area (yellow). Assumptions about fabrication, coating, and seeds of the scaffold. In the currentstudy, processing a scaffold relates to setting the model parameters. There are several existing techniques available for fabrication and modification of the scaffold, for example the poly (D,L-lactide-co-glycolide) (PLG) with a gas-foaming/particulate-leaching process [14,24]. To prepare the coating and seeds of the scaffold, hydroxylapatite (HA) nanoparticles (Sigma-Aldrich) were suspended in phosphate-buffered saline and sonicated for 1 min to dissociate the aggregates [14,24]. The nanoparticles were then complexed with several ECM components (forming HA/ECMs), including collagen I (BD Biosciences, San Jose, CA), fibronectin (Sigma-Aldrich, St. Louis, MO), and laminin I (Trevigen, Gaithersburg, MD) [22]. Lentivirus was used to encode the neurotrophic factors NT-3 (HA/LV-NT-3) and BDNF (HA/LV-BDNF) [24] and the ChABC gene (HA/LV-ChABC) [21,23,24]. The HA/ECMs and HA/virus complexes were then incubated for 10 min at 4˚C and deposited onto the scaffold surface using Stripper pipette tips (Mid-Atlantic Diagnostics, Mount Laurel, NJ). The HA/ECMs and HA/LV-ChABC were coated over the entire surface, whereas HA/LV-NT-3 and HA/LV-BDNF were only deposited on theoff-ramp surface of the bridge. Although HA could firm the coating and/or seeds on the PLG material surface [24], it probably stimulates bone formation [34]; therefore, it would be best to avoid the overuse of HA. In addition, the prepared scaffold should be placed on dry ice until implantation. This scaffold could then beused to replace previous scaffolds used in a rat model for SCI repair [10,11,16], to bridgea gap of approximately 3 mm. The scaffold was considered to perform the following functions:the factors seeded on the bridge surface were localized and sustained; the on-ramp slope introduced the rostral regenerative axons onto the scaffold; the secreted ChABC cleared the growth pathways; the ECM components kept the axons on target; and the gradients of the secreted NT-3 and BDNF guided the growing rostral axons along the correct path while blocking the caudal axons.

Scaffold surface equation and constraints for growth cones
Geometrically, the scaffold is an ellipsoid that rotates about its z-axis ( Fig 1B). The surface equation reads as follows: Here, x, y, and z are the coordinates of a point on the surfaceand R is the distance from that point to the z-axis.
When an axon grows along the scaffold surface, the growth cone can be considered as a particle moving on the surface. The ECM molecules coated on the scaffold are assumed to strongly adhere to the membrane proteins of the growth cones, tethering them to the scaffold. In other words, the movements of the growth cones are constrained by Eq (1). Taking the time derivative of both sides of Eq (1), the constraint conditions of the growth cone velocity can then be obtained as follows: where V x , V y , and V z are the velocity components of a growth cone in the x-, y-, and z-directions, respectively, and V R is the velocity component along the rotational radius R. In Eq (2), V z is the longitudinal velocity of the growth cone, which defines the drawing speed or growth rate of the axons. When V z >0, the axon is elongating. Otherwise it is retracting. V R is the lateral velocity of the growth cone. When V z >0 and z<0, V R >0 and the growth cone progresses from the on-ramp of the scaffold. Conversely, when V z >0 and z>0, V R <0 and the growth cone advances to the off-ramp. In the next two subsections, we describe the drawing force and locomotive guidance of the growth cone along the scaffold.

Equations for axonal growth
Many pieces of evidence have shown that axonal growth cones move because of chemotactic processes, biased toward (attractive chemotaxis) or away from (repulsive chemotaxis) the chemical source [35,36,37]. The chemicals that attract and repel regenerative axons are found on and around the multifunctional scaffold (this is hereafter assumed to be the case in the implanted status), and will be classified in the next subsection. The chemotactic force, which defines the attractive or repulsive action on a growth cone, is proportional to the gradient of the diffusible molecules released from the chemical source [29,32,37,38,39,40,41]. The growth cone can be regarded as a particle whose persistent velocity determines the growth rate of an axon. As axonal growth is particularly slow (approximately 0.01 μms −1 ) [35,36,37], the acceleration or inertial force can be neglected [24]. Stochastic factors can also be neglected because the chemotactic movement is highly consistent [29]. Therefore, the velocity of the growth cone is directly proportional to the chemotactic force: where k and N A are the number and total number of regenerative axons/growth cones, respec- where i, j, and kare the unit vectors in the Cartesian coordinate system) is the position of the k-th growth cone at time t. μ is the dynamic viscosity coefficient.
is the result of the chemotactic forces from all types of chemotacticrelated molecules (CRMs) acting on the k-th growth cone at t. The CRMs will be classified into three types in the next subsection, where i is the number of types of CRMs, and subscript idenotes the i-th type of CRMs (CRMs-i). The chemotactic force [Eq (5)] is defined as a dimensionless vector (p i ) with a proportionality constant (λ i ). In this equation, ρ i is the concentration of CRMs-i at r A k and t; rρ i is the gradient of ρ i (where r = i@/@x+j@/@y+k@/@z is the Hamiltonian operator); jDr A k j ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi Dx 2 þ Dy 2 þ Dy 2 p is the scalar difference of r A k across the width of the k-th growth cone; and ρ ∑ is the sum of ρ i over all types of CRMs. Note that the scalar of p i [Eq (5)] simplifies to Δρ i /ρ ∑ , which expresses the relative difference of ρ i , with Δρ i being the absolute difference of ρ i across the distance jDr A k j. In practice, the average width of the growth cone can be considered as jDr A k j~10μm [29,38]. Therefore, the gradient and relative difference of the chemotactic concentration are mathematically linked through p i . The corresponding proportionality constant λ i then acquires a clear physical meaning of force per unit length. This model considerably reduces the growth-rate distortion of the growth cone when close to the target cells. Note that in the one-dimensional single-component problem [29,38,39], p i reduces to p = Δρ/ρ, which defines the concentration gradient in some biophysical areas.
Finally, from Eqs (2), (4) and (5), the drawing speed V z of the k-th axon can be expressed as follows: where Δz is the cone's width along the z-axis. The other symbols have been defined in Eqs (4) and (5).

Numerical methods
Eqs (1)-(10) comprise the mathematical model of regenerative axons growing along the spherical multifunctional scaffold where their status is as shown in Fig 1A or Fig 2A-2C. The source terms in Eqs (9) and (10) are nonlinearly coupled with Eq (8) and Eqs (4)-(7) through point r A k that tracks the growth cone. Therefore, Eqs (4)-(10) comprise a set of coupled nonlinear differential equations that can only be solved numerically.
In the numerical simulation, we changed the lateral size of the scaffold (i.e., we changed the slope of the on-ramp/off-ramp) and the seeding densities of the HA/ECM components HA/ LV-ChABC, HA/LV-NT-3, and HA/LV-BDNF. The simulation recorded the growth rates of the regenerating axons, and tested whether the regenerating axons could grow across the scaffold and finally connect with their targets.

Consideration of the parameter values
Based on data from in vitro experiments [35,36,37], we can estimate the order of magnitude of each parameter. The estimates are as follows: growth cone width jDr A k j is approximately 10-20μm; axonal growth rate is approximately 0.01μms -1 ;and the diffusion coefficient D 1 and dissociation constant K d of CRMs-1 (e.g. NGFs) are approximately 10-50μm 2 s -1 and 1nM, respectively. In addition, the range of concentrations for NGFs is~(0.01-10)K d , and the minimum relative concentration difference to which the growth cone can respond [35,36,37,38,39], |Δρ i /ρ i |, is~1%. However, the values of many parameters cannot be recognized, including the point-source release rate σ i , the attenuation coefficients k −i , the force proportionality constant λ i , and the viscosity coefficient μ. Given that the diffusion velocity k À 1 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi of CRMs-1 should exceed the constant velocity of the growth cone, the absolute concentration should satisfy ρ 1 �K d and the relative concentration difference should satisfy ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi D 1 =k À 1 p � 1% at the diffusion radius ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi . Under these conditions, the attenuation coefficient was estimated as k −1 �2.0×10 −6 . At the point source, we should have ρ 1 = eK d +σ 1 /k −1 �10 K d so that we can estimate the point-source release rate as σ 1 radius was then estimated as r max = 5.6 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi . Finally, in the concentration field based on the above data, the axonal growth rate was assumed to be 0.0025-0.05μms -1 . After many numerical trials, the number of λ 1 /μ was approximated as 1. The parameter values of CRMs-2 and CRMs-3 were set based on those of CRMs-1. All parameter values are listed in Table 1. To save computational resources, we employed coarse-grained processing, i.e., making up an equivalent number of seeds (point sources) for calculation that was smaller than that of seeding on the scaffold during fabrication. Because the fabrication process controls the densities of HA/ECMs/LV-ChABC, HA/LV-NT-3, and HA/LV-BDNF seeded on the scaffold, the densities of the point sources of CRMs-1, CRMs-2, and CRMs-3 can be quantified as , respectively, where g 0 1 , g 0 2 , and g 0 3 are the normal point-source densities of CRMs-1, CRMs-2, and CRMs-3 (in μms -2 ) respectively, and η 1 , η 2 , and η 3 are their respective controllable scale factors. The areas of the scaffold off-ramp (caudal side) and the ventral portion of a growth cone are denoted by A C and A cone , respectively. Thus, the total number of CRMs-1 point sources in the off-ramp area is γ 1 A C , i.e., N T = γ 1 A C in Eq (8). Meanwhile, the total numbers of CRMs-2 and CRMs-3 connected to a single growth cone are γ 2 A cone and γ 3 A cone , respectively. In Eqs (9) and (10), N A is the total number of axons severed in an injury event. In a complete transection rat model, N A might exceed 10 4 (estimated from the cross-sectional area ratio of an axon to the spinal cord). Because the exact value N A was unavailable, we ranged N A from 300 to 12000 in the simulations. Correspondingly, N T and σ 1 in Eq (8) were replaced by their equivalent values N eqT and σ eq1 , respectively. Considering N eqT = N A and applying the mass conservation principle, we obtained σ eq1 = γ 1 A C σ 1 /N A , where A C /N A = A C /N eqT represents the area occupied by each equivalent point source. Here, the point sources were assumed to be evenly distributed on the off-ramp of the scaffold. Similarly, σ 2 and σ 3 in Eqs (9) and (10) were replaced by their equivalent values σ eq2 = γ 2 A cone σ 2 and σ eq3 = γ 3 A cone σ 3 , respectively, where A cone represents the area on the scaffold surface covered by one moving growth cone at that moment and γ 2 A cone and γ 3 A cone represent the original numbers of CRMs-2 and CRMs-3 point sources on one A cone , respectively. After an equivalent transformation, each A cone includes two point sources: a CRMs-2 point source with release rate σ eq2 and a CRMs-3 point source with release rate σ eq3 . Both releases activate at the location of the cone and deactivate when the cone departs. This activation-deactivation process of the point sources is modeled by the Dirac delta function δ(�) in Eqs (8)- (10).
In the simulations, we confined the architecture and mathematics of the model to a cubical compartment with side length L = 5000 μm and imposed absolute boundary conditions. The D3Q15 mode [32,47] with a 64 ×64 ×64 lattice was applied to LBM simulations. MATLAB 7.11.0 (The MathWorks, Inc.) was employed to programthesimulations and run them on ThinkServer TD350 (Lenovo Group Ltd). All parameters used in the simulations are summarized in Table 1.

Results
Before describing the simulation results, we provide the following definitions. A severed axon is considered to have successfully regenerated when it has regrown along the scaffold surface from the on-ramp to the off-ramp within 2 weeks of treatment. The success rate of axonal regeneration is the number ratio of successfully regenerated axons to all severed axons in the injury event (calculated by the current mathematical model). The growth rate of the from the lower part (the rostral/on-ramp) of each scaffold represent the regenerated axons (numbering 0 .78 N A , 0.7 N A , and 0.58 N A in a, b, and c, respectively, corresponding to axonal regeneration success rates of 0.78, 0.7, and 0.58, respectively). d)-g) Results related to a typical regenerating axon taken from the calculation of scaffold No. 1 (qualitatively similar results are obtained for the other two scaffolds). At the beginning, when the axon is at the on-ramp, its growth is unstable and some retraction is observed (d and e) because the concentrations of the positive factors CRMs-1 (f) and CRMs-3 are low (g; green), while those of the inhibitoryfactors are high (g; red). This situation is reversed as the axon approached the off-ramp.
. The growth rate should exceed 0.0025 μms -1 (1.5 mm per week) but should not exceed the physiological limit (~0.05 μms -1 ). The effective growth rate of the axons is the product of their success and growth rates, i.e., the longitudinal extension velocity averaged over all severed axons in the injury and repair process. It is given by

Influence of the number of severed axons on regeneration
The number of severed axons N A usually depends on the damage event. As previously mentioned, the maximum number of severed axons might be in the order of 10 −4 . In the following computer simulations, we varied N A from 300 to 12000 and observed the consequent changes and determined the most efficient size for axonal regeneration. For this purpose, we investigated how the slope of the on-ramp/off-ramp of the scaffold affects axonal regeneration. In this section, the density scale factors of the CRMs-1, CRMs-2, and CRMs-3 point sources on the scaffolds were considered as η 1 = η 2 = η 3 = 1. In other words, the seed densities were constant and set to the same value (specifically, γ 1 = γ 2 = γ 3 = 0.1 μm −2 ; see Table 1 and the subsection "Equivalent number of seeds on the scaffold and the point-source release rate" for details). All other parameter values are listed in Table 1. To clarify the scope and level of the model and the developed method, we first discuss the special case of N A = 300. Fig 2 demonstrates the growth path and velocity of the regenerative axons and the CRM concentrations around the growth cone. By increasing N A from 300 to 12000 at irregular intervals, we obtained a series of results similar to those demonstrated in Fig 2. These results are shown in Fig 3, which reveals how the number of severed axons influences axonal regeneration. Fig 3A shows the effect of the number of severed axons N A on the growth rate of the regenerating axons on each scaffold. Each point in the curves indicates the average growth rate of the surviving axons that grew along the scaffold from the on-ramp to the off-ramp. Axons that died or failed halfway were excluded. The error bars indicate the standard deviations of the averages. N A does not significantly affect the growth rate provided that the axons can successfully regenerate. Moreover, the growth rate is proportional to the lateral size of the scaffold.
In Fig 3B, the success rate of axonal regeneration decreases with increasing N A . When N A was moderately large (<9000), the slim scaffold achieved a high success rate. When N A = 300, the success rate was approximately 35% higher on scaffold No. 1 than on scaffold No. 3; however, when N A = 12000, the lateral size effect of the scaffold disappeared and the success rates on all three scaffolds declined to <30%. As the number of severed axons increased from 300 to 12000 on scaffold No. 1, the success rate was reduced by a factor of 3. From Fig 3A and 3B, we infer that a small on-ramp slope (i.e., a slim scaffold) increases the success rate of axonal regeneration up to a certain number of severed axons. When too many regenerative axons pass over the narrow bridge, the congestion lowers the growth rate. If success rate is more important than growth rate, the slim scaffold would be the first choice. Fig 3C shows the effective growth rate of the axons (the product of the growth rate in Fig  3A and the success rate in Fig 3B as N A increases). Each point in the curves represents the growth rate of all regenerated axons, averaged over all severed axons, in the injury and repair process. The effective growth rate is a comprehensive index of axial sprouting and growth of axons after an SCI. Again, the slim scaffold is advantageous for axonal regeneration when N A �6000.
In Fig 3D, the concentration of CRMs-1 sensed by the regenerating axons was relatively robust to both N A and the lateral size r a of the scaffold. However, the concentration ratio of CRMs-3 to CRMs-2 sensed by the regenerating axons heavily depended on both N A and r a (Fig 3E). Note that N A and r a are also highly correlated with the success rate of axonal regeneration ( Fig 3B) and the effective axonal growth rate (Fig 3C). This suggests that the concentration ratio of CRMs-3 to CRMs-2 deserves more attention than it currently receives in SCI treatments. Fig 3F averages the longitudinal coordinates in Fig 3B over   Numerical characterization of a new multifunctional scaffold for SCI

Influence of CRMs-1 point-source density on regeneration
This analysis was performed on scaffold No. 1 (r a ×r b = 0.15×0.3). The prototype of the CRMs-1 point source is HA/LV-NT-3 and HA/LV-BDNF seeded on the off-ramp (Fig 1B,yellow) of the scaffold during fabrication (see subsection "Coating and seeds for the scaffold"). In each test, the seeding density of CRMs-1was g 1 ¼ Z 1 g 0 1 (g 0 1 = 0.1 μm −2 ; see Table 1), where η 1 was increased from 0 to 100 at irregular intervals. Alternatively, the CRMs-2 point source includes various inhibitory components, for example MAG induced by debris from the myelin sheath resulting from when axons were damaged, and CSPGs remaining on the scaffold surface, which cannot be easily controlled. Here, the density of the CRMs-2 point source was graduated through six levels: g 2 ¼ Z 2 g 0 2 (g 0 2 = g 0 1 = 0.1 μm −2 , with η 2 = 0, 1, 5, 10, 50, and 100). The density of the CRMs-3 point source was then expressed as g 3 ¼ Z 3 g 0 3 . CRMs-3 represents the HA/ECM components (collagen I, fibronectin, and laminin I) coated on the whole surface of the scaffold during fabrication. Therefore, we set η 3 = 1 (or g 3 ¼ g 0 3 = 0.1 μm −2 , a constant). The number of severed axons can be selected from the feasible range (300-12000). Considering the available computer resources, we maintained N A to be constant at 1200. Under these conditions, we performed the calculation and obtained Fig 4. Fig 4A shows the growth rates of the regenerating axons as functions of density (η 1 ) of the CRMs-1 point source for different densities (η 2 ) of the CRMs-2 point source. Axonal growth is driven by the concentration gradient of CRMs-1 (rρ 1 ) and is supported and inhibited by rρ 3 and rρ 2 , respectively. Note that ρ 2 and ρ 3 are coupled to ρ 1 through Eqs (8)- (10). Therefore, the growth rates of the regenerating axons are not simply proportional to η 1 but change nonlinearly with increasing η 1 . Moreover, their extreme values depend on η 2 . Note that the growth rate was always lowest when η 2 = 0, i.e., when no inhibitors were present in the injured microenvironment. This suggests that some remaining inhibitors (e.g., 0<η 2 �1) are required to increase the growth rate. However, when η 1 >15 and η 2 >5, the growth rates were quite similar, suggesting that once the inhibitor level has exceeded a certain threshold, simply increasing the promoter levels is ineffective for increasing the growth rate. Fig 4B shows the success rate of axon regeneration with increasing density (η 1 ) of the CRMs-1 point source. First, when η 2 = 0, i.e., when no inhibitor was present in the injured microenvironment, the success rate was 100% for all η 1 , except when η 1 = 0 (i.e., when the success rate was 0). Next, for any fixed η 1 >0, increasing η 2 decreased the success rate. Finally, when η 1 was sufficiently high, the success rate of axonal regeneration exceeded 92%, regardless of η 2 . This suggests that reducing the inhibitor while increasing the promoter improves the success rate of axonal regeneration. The latter may be more important than the former because even if η 2 reached 0, η 1 >0 would be required for axonal growth. However, as previously mentioned, a high success rate implies overcrowding of the growing axons on the scaffold, which eventually causes congestion and lowers the axonal growth rate. If the growth rate is excessively low, the axons might abort the growth process (in practice, growth cessation is due to unexpected causes) or an opportunistic time window for subsequent treatments may be missed. Therefore, when predicting axonal regeneration, the success rate must be balanced against the growth rate of the regenerating axons. Fig 4C shows the effective growth rate of the axons. As previously mentioned, this index represents the sprouting and growth of the injured axons. Additionally, it is the product of the growth rates in Fig 4A and the success rates in Fig 4B. The effective growth rates were higher for η 1 >5 and η 2 = 1 than for η 2 = 0. If η 1 >50 could be achieved in practice, three treatments (η 1 >5and η 2 = 1, η 1 >30 and η 2 = 5, and η 1 >50 and η 2 = 10) are preferred over the η 2 = 0 treatment, and increasing η 1 is much easier than achieving η 2 = 0. Fig 4D shows the average success rates of axonal regeneration over all η 1 (0-100) for each η 2 . Regardless of η 1 , the success rates on scaffold No. 1 decreased with increasing η 2 and the marginal effect was large for small η 2 . Again, this result highlights that reducing the inhibitors in the microenvironment is important for axonal regeneration. Note that when η 2 = 0, the average success rate should be less than one (below 100%) because the success rate at η 1 = 0 was zero. Thus, its logarithm (horizontal axis in Fig 4B) could not be defined. However, owing to non-zero η 1 a unity success rate was achieved (Fig 4B).

Influence of CRMs-3 on axonal regeneration
Finally, we numerically test the hypothesis that an over-eutrophic scaffold surface harms axonal regeneration.
The analysis was performed on scaffold No. 1 (r a ×r b = 0.15×0.3), and assumed a constant number of severed axons (N A = 1200). The prototype of the CRMs-3 point source comprises HA/ECM components (collagen I, fibronectin, and laminin I) coated on the whole surface of the scaffold, as described previously. The coating density of CRMs-3, expressed as g 3 ¼ Z 3 g 0 3 with g 0 3 = 0.1 μm −2 (see Table 1), can be set during fabrication. In each test, density η 3 was increased from 0 to 100 at irregular intervals. The CRMs-1 and CMRs-2 point-source densities were set as g 1 ¼ Z 1 g 0 1 (with g 0 1 = 0.1 μm −2 and η 1 = 1) and g 2 ¼ Z 2 g 0 2 (with g 0 2 = 0.1 μm −2 and η 2 = 0, 1, 5, 10, 50, and 100), respectively. The calculation results under these conditions are presented in Fig 5. As shown in Fig 5A, the growth rates of the regenerating axons decreased with increasing η 3 of the CRMs-3 point-source density, regardless of η 2 . Over a wide range of lower η 3 values, the growth rates remained stable for all η 2 , except for η 2 = 0. However, beyond a threshold η 3 value, which increased with increasing η 2 , growth sharply declined. This supports the hypothesis that an over-eutrophic scaffold surface impedes axonal growth, and that certain inhibitors might naturalize the over-eutrophic effect and weaken or delay the harm. Fig 5B shows the success rates of axonal regeneration varying as functions of η 3 for each η 2 . When η 2 = 0, i.e., when no inhibitor was present in the injured microenvironment, the success rate remained at 100% till the point when η 3 = 60, and then dropped sharply to 16.67% at η 3 = 70. Subsequently, the success rate fell to zero. At non-zero η 2 , the success rates first slowly ascended with increasing η 3 , then rose sharply before dropping to a low level, and finally reached zero. When η 2 was high, the success rate curve began at a low level and ascended very slowly, forming a low plateau. For η 2 �50, the success rate peaked only when η 3 exceeded 100. This suggests that the use of an over-eutrophic scaffold surface to improve the success rate of axonal regeneration is inefficient and risks destabilization of the regeneration process. Fig 5C plots the effective growth rate of the axons as functions of η 3 for each η 2 . Applying an over-eutrophic scaffold surface to promote the sprouting and growth of axons after an SCI might be inefficient or even counterproductive. Fig 5D shows the average success rates of axonal regeneration over η 3 = 0−100 for each η 2 . Regardless of η 3 , the success rates decreased with increasing η 2 on identical scaffolds and the marginal effect was greater at low η 2 than at high η 2 . As discussed earlier, this result highlights the fact that reducing inhibitors in the microenvironment is important for axonal regeneration.

Discussion
This study specifically aimed to create a microenvironment for axonal regrowth after SCI by mathematically changing the location of seeding on the scaffold, the density of the cells/factors seeded, and the size and shape of the scaffold. Therefore, we assumed a solid, spherical, multifunctional, biomaterial scaffold that bridges the rostral and caudal stumps of a completely transected spinal cord in a rat model for the calculations.
We assumed the body of the scaffold was made of PLG, the whole surface was coated with HA/ECMs and HA/LV-ChABC, and its off-ramp at the caudal area was additionally seeded with HA/LV-NT-3 and HA/LV-BDNF. These factors perform several functions. First, the onramp slope at the rostral area of the scaffold steers the growth cones of the regenerative axons onto the scaffold smoothly. Second, the factors seeded on the scaffold surface can be localized and sustained over a reasonably long period of time. The HA/ECM components adhere to the growth cones on the scaffold, supporting regenerative axons. HA/LV-ChABC secretes ChABC, which degrades CSPGs, inhibitory components from the glial scarsurrounding the injured tissue. The NT-3 and BDNF molecules released from HA/LV-NT-3 and HA/ LV-BDNF at the off-ramp area can diffuse to the on-ramp area; this gradient can then guide the axonsat the rostral stump to grow along the outer surface of the scaffold and enter the caudal stump area, while also preventing the axons regenerating at the caudal stump from growing toward the rostral stump until they connect with axons emerging from the on-ramp area. That is, the scaffold forms a one-way bridge from the rostral to the caudal side, with a gentle slope of entry. From a mathematical point of view, the profile curve of the scaffold results in a small and continuous tangent slope, and the resulting concentration of all factors (promoters and inhibitors) on and around the scaffold, and varying along the span of the scaffold from one side to another, causes a monotonic increase. Provided that the scaffold is implanted upside down, the gradient of the resulting concentration or the direction of growth of axons will be reversed correspondingly, inferring from previous observations [18] that axons grow toward a Numerical characterization of a new multifunctional scaffold for SCI site where additional cells have been injected, so either side of the bridge could be used as an injection site.Note that no matter which side of the bridge is used as the entry point, both motor and sensory neurons axons can cross the bridge because axons of both motor and sensory origins have previously been found within a tunneled scaffold [14].
Several scaffolds for SCI have previously been described [3,4,5,6,7], such as cell grafts [8,9], tubes or conduits [10], and cylinders with tunnels/linear pores or filled with fibers, either randomly or in alignment. Of these, cell grafts comprise natural soft tissues, with conduits that are empty or contain soft matrices and/or cells, and with a wide diameter, making it easy for regenerative axons to enter them, with no need to consider an entry slope. However, axons which entered these grafts were often trapped in them and rarely re-entered the host tissue. It was commonly believed that fine physical guidance was required within the cell grafts and conduits. Therefore, cylinders with tracks (a rolled-up nanofiber sheet or film) [10,11] or tunnels (via modeling) [14,19], combined with seeding of cells/factors [12,13,14,15,16,18,19], emerged and mostly replaced the use of cell grafts and conduits. However, not only did regenerative axons still get trapped in the guideways [10,11,12,13,14,15,16], but congestion was also observed at the entries [19]. This trapping, according to the current study, is caused by the over-eutrophication in and/or on the scaffold, and the lack of additional chemotactic factors close to either side of the scaffold. Even if these factors are uniformly distributed throughout whole scaffold during fabrication, because the factors diffuse easily from the two ends of the scaffold, they reach the middle of the scaffold after a period of time, which attracts regenerating axons from both ends to the middle. Therefore, the scaffold becomes a two-way channel, leading the axons to grow toward each other, probably resulting in their intersection. However, there is little direct evidence to show that such intersected axons can form synaptic connections. In this situation,the factors on the scaffold are more likely to be a barrier against the regenerating axons. Fortunately, the axons can sometimes break through this barrier [19], as was simulated in this study, and re-enter the host tissue. This is typically due to additional chemotactic factors being set close to either the caudal [19] or rostral stumps [18]. Alternatively, factors (known or unknown) contained in and/or on the scaffold might form (intentionally or unintentionally) a consistent gradient along the span of the bridge; there may also be a small probability that some unknown factors from host tissues form a beneficial gradient for axonal regeneration. Based on the typical situation [19], we inferred that, according to the principle of chemotaxis of axons, only if the peak of the barrier is much lower than the summit formed by the chemotactic factors at either end, the axons can break through the barrier.
However, another barrier that exists for scaffolds with tracks or tunnels is the entry obstruction that lowers the number of axons entering into the spaces or pores of the scaffoldbecause of the absence of a uniform entry slope. For tunneled scaffolds [14,19], the slope of entry is small or zero only for those axons whose growth cones face to the pores, where the axonal incidence angle α = 0 (i.e., the slope = tan(0) = 0); whereas for the other axons at the entry, the slope is abrupt or infinite, due to α = ±π/2 (i.e., slope = tan(±π/2) = ±1). This prevents the axons from regenerating straight ahead and reduces the number of axons entering the tunnels. Pore size is another parameter that affects the ability of axons to pass through a scaffold, and while investigating optimal pore size is difficult, the present model can be modified to study this. Once the number of axons entering a scaffold is small, the number of axons exiting will be much lower, according to previous observations [19]. The present spherical scaffold described by us is not only simple but also, at least theoretically, lacks the shortcomings inherent in other porous scaffolds.
The proposed scaffold could be used to replace previous scaffolds used in rat models for SCI repair [10,11,16], to bridgea gap of approximately 3 mm. It is worth noting that the present scaffold might offer additional possibilities. The complexesseeded on the caudal area that express NT-3 and BDNF, as well as being diffusible chemoattractants for attracting axonal growth cones, have enhanced oligodendrocyte survival and axon myelination [24,51,52,53]. NT-3 and BDNF may also enhance a number of other processes. For instance, BDNF has been associated with a reduction in the inflammatory response, including a reduction in astrocyte numbers [54], which further aids axon regeneration. Whilefabrication of our scaffold might be difficult, and the use of hydroxylapatite (HA) is probably not the best choice due to the risk ofMilwaukee shoulder [34], we do provide a prospective direction for SCI repair.
Our mathematical model embodies the geometry, chemistry, and physics of the system under investigation. The shape and size of the scaffold provide the geometrical boundary conditions that constrain the growth cone's movement or the regenerative axon growth. While the inclusion of all chemical factors in a single model is difficult, a coarse-grained method can be used to obtain a balance between the reduction in the number of factors and the retention of the chemical properties. That is, the CRMs on and around the scaffold (assuming it has been implanted) were classified into three types with different chemical properties: the CRMs-1 group comprised chemoattractants of axonal growth (NT-3/BDNF secreted by seeded HA/ LV-NT-3 and HA/LV-BDNF at the off-ramp area); the CRMs-2 group comprised chemorepellents of axonal growth (compounds produced by the injured tissue, such as Nogo-60, MAG, and OMG, and the remnant CSPGs that are not neutralized by ChABC released from seeded HA/LV-ChABC); and the CRMs-3 group comprised molecules released from the coated HA/ ECM components, which support axonal growth. Among these groups, CRMs-1 plays the leading role in axonal regeneration, whereas CRMs-2 and CRMs-3 provide a balanced and coordinated effect, and interact with CRMs-1 (for which, CRMs-2/3 was formulated as a function of CRMs-1). The physical and biophysical aspects of the model are the Fickian diffusions and reactions of the CRMs and the chemotaxis of the axonal growth cone motility. Under CRMs-1 and CRMs-3 gradients, the regenerative axons elongate toward the target cells, whereas under the CRMs-2 gradient, they retract.A similar chemotaxis pattern has been validated via experiments [35,36,37], and mathematicallymodeledfor studying axonal growth in neural development [28,29,30,31]. Note that axonal growth/regrowth follows the same processes observed experimentally both in development and injury, in vivo and in vitro, because in both cases the experiments always initially cause damage to the nerve cells, for instance during surgery or separation. The difference between them, however, is in the degree of damage, and in the age of the experimental subjects. Therefore, a theoretical model for axonal growth during development can be modified for studying SCI. In fact, the present model is mathematically similar to those described in the literature [28,29,30,31].
By applying our mathematical model to the theorized scaffold, we numerically studied the influence of the number of severed axons; the slope of the on-ramp of the scaffold (which is related to the lateral size of the scaffold); and the concentrations and gradients of the CRMs (which are related by their seeding densities to the survival and growth of the regenerative axons).
Axonal regeneration was evaluated based on the growth and success rates of the regenerated axons. A severed axon has successfully regenerated when it has regrown along the scaffold surface from the on-ramp to the off-ramp within 2 weeks of the treatment. The success rate defines the number ratio of the successfully regenerated axons to all axons severed in an injury event. The growth rate is the average longitudinal extension velocity of the successfully regenerated axons.
It should be noted, however, that our model is highly idealized, for example, the basis data we used for the calculation of growth cone velocity were not derived from trauma tissue, which might be spatially different in physics and chemistry, and not be reflected by CRMs-2. The intersections and/or tangles between regenerating axons and the energy expenditure for axonal growth were not modeled in the current study. The cut to generate a spinal cord injury should not be too wide (>1 cm) for using this model to predict SCI because the effective diffusion distance of CRMs-1 might be limited to approximately 1 cm by Fick's first law, on which this model is based. Only an in vivo verification will show whether these equations will hold.

Conclusions
A solid, spherical, multifunctional, biomaterial scaffold is assumed to bridge the rostral and caudal stumps of the spinal cord in a completely transected rat model, thereby promoting the entry of regenerative axons from the rostral stump into the caudal stump tissue at the opposite side of the scaffold.
Three scaffold shapes (slim, round, and stocky) were investigated in our simulations. Among them, the slim shape benefited axonal regeneration the most by presenting a small slope at the on-ramp area. However, if the success rate becomes too high, numerous regenerative axons crowd into a narrow area, causing congestion and resulting in a reduced growth rate. The stocky scaffold induced the opposite effect, and the round scaffold induced intermediate effects. When success rate is more important than growth rate, the slim scaffold should be the first choice.
The number of severed axons in an injury event (between 300 and 12000) does not significantly affect the growth rate of the regenerated axons, but does influence the success rate of axonal regeneration (particularly, the success rate decreases with increasing number of severed axons).
Among the three types of chemical treatments, raising the CRMs-1 (NT-3 and BDNF) level while reducing the CRMs-2 level (CSPGs and other chemorepellents) benefited the success and growth rates of axonal regeneration the most. Physically, the CRMs-1 level was increased by increasing the seeding density of HA/LV-NT-3 and HA/LV-BDNF on the off-ramp of the scaffold, whereas the CRMs-2 level was reduced by increasing the seeding density of HA/ LV-ChABC over the entire scaffold surface. However, raising the CRMs-3 (ECM components) level by increasing the density of HA/ECM components over the entire scaffold surface may create an over-eutrophic surface that harms axonal regeneration.
The theoretical predictions made in this study need to be experimentally validated in the future. In principle, the current tool can be easily modified for predictions regarding scaffolds with other architectures.