How a reaction-diffusion signal can control spinal cord regeneration in axolotls: A modeling study

Summary Axolotls are uniquely able to completely regenerate the spinal cord after amputation. The underlying governing mechanisms of this regenerative response have not yet been fully elucidated. We previously found that spinal cord regeneration is mainly driven by cell-cycle acceleration of ependymal cells, recruited by a hypothetical signal propagating from the injury. However, the nature of the signal and its propagation remain unknown. In this theoretical study, we investigated whether the regeneration-inducing signal can follow a reaction-diffusion process. We developed a computational model, validated it with experimental data, and showed that the signal dynamics can be understood in terms of reaction-diffusion mechanism. By developing a theory of the regenerating outgrowth in the limit of fast reaction-diffusion, we demonstrate that control of regenerative response solely relies on cell-to-signal sensitivity and the signal reaction-diffusion characteristic length. This study lays foundations for further identification of the signal controlling regeneration of the spinal cord.


Validation of the modelling of the spatial motion of the ependymal cells
In this section, we aim to validate the cellular scale of our model (described in the STAR Methods section of the main text) by using the 1D model developed in 1 for spinal cord regeneration.Indeed, the model of 1 bears analogies with our own model as cell division follows the same rules as described in the STAR Methods section of the main text and cells are supposed to behave as rigid spheres.However, cell motion is treated differently in 1 as upon division, one daughter cell occupies the same position of the mother cell and the other daughter cell is placed posteriorly to the mother cell, moving all the posterior cells a quantity corresponding to the size of a cell in the posterior direction.This corresponds exactly to the dynamics of rigid spheres in a 1D setting, but does not extend easily to 2D, justifying the need for the more complex treatment of cell motion that we proposed (see STAR Methods section of the main text).The goal of this section is therefore to validate our treatment of cell motion by comparing it to the exact motion of rigid spheres provided in 1 .
To this aim, we also need to adapt our signalling process, as in 1 the ependymal cells are supposed to be recruited by a signal spreading with constant velocity from the amputation plane, contrary to our more realistic model of a signal evolving via a reaction diffusion equation (see STAR Methods section of the main text).Therefore in this section, we consider the same signal as in 1 , i.e a signal moving with constant velocity from the amputation plane, during τ = 85 ± 12 hours and recruiting ependymal cells within λ = 828 ± 30 µm from the amputation plane.
In Fig. S1 C, we show the tissue outgrowth (distance from the amputation plane to the the most posterior ependymal cell) obtained with our modified model in 1D (red curve) compared to the outgrowth of with the same set of parameters as in 1 .Our model prediction is in good agreement with the previous results 1 , validating the cellular scale of our model in 1D.

Validation of the numerical scheme for the reaction-diffusion process of the signal
In this section, we aim to validate the numerical scheme of the reaction-diffusion process governing the signaling scale of our model (see STAR Methods section of the main text).To this aim, we consider a simplified model where the tissue is fixed (no cell division), and 1D.Under these conditions, the equation for the signal simplifies to a reaction-diffusion equation on a fixed 1D domain with Dirichlet boundary conditions on one side and no-flux boundary conditions on the other side.This enables to derive an analytical solution adapting the method used in ceccarelli2021 to our setting, and compare our numerical simulations with the exact solution obtained.Derivation of the analytical solution of the reaction-diffusion model in finite domain.
Here, we derive an analytical solution for a reaction-diffusion equation in a finite and fixed 1D domain, with non-homogeneous Dirichlet conditions on one side and Neumann boundary conditions (null flux) on the other side.We adapt the derivation of ceccarelli2021 , where the authors considered various but different boundary conditions.
We consider a signal spreading in a non-growing 1D tissue of length L and we call C 1 (z, t) its concentration, with z ∈ [0, L] and t ≥ 0. In order to adapt the derivation of ceccarelli2021 , we consider a new spatial variable x = L − z and a 1D tissue of length L starting at x = 0. We assume that the signal concentration is given at x = 0, it diffuses to the tip of the tissue with a diffusion constant D and degrades linearly at a rate k.We assume no-flux boundary conditions at the tip of the tissue in x = L.At t = 0 there is no signal in the tissue.The evolution in space and time of the signal distribution C 1 is described by the following reaction-diffusion equation We have the following initial and boundary conditions: • no signal at initial time: • fixed signal concentration at x = 0: where C 0 is constant, • Neumann boundary conditions at the tip of the tissue x = L: We introduce the dimensionless variables where λ = D k is the characteristic reaction-diffusion length, and we rewrite Eq. (S.1).We define the quantity ℓ = L λ , which is now the only model parameter.We rewrite the concentration as and we obtain the following equation with the following initial and boundary conditions: • no morphogen at initial time: • fixed signal concentration at ε = 0: where C 0 is constant, • Neumann boundary conditions at the tip of the tissue ε = ℓ: To solve this equation we introduce an auxiliary function C 2 and we define C as follows: In order to rewrite Eq. (S.3) in terms of C 2 , we compute the spatial and time derivatives of C as follows: ∂ε 2 and ∂C ∂τ computed in Eq. (S.4) and Eq.(S.5), respectively, into Eq.(S.3), leads to the following equation with the following conditions To solve Eq. (S.6) we have to redefine C 2 using the auxiliary functions C 2 (ε, τ ), f (τ ) and g(ε).We define while f (τ ) and g(ε) will be defined explicitly later on.The time derivative of C 2 (ε, τ ) is and the spatial derivative is Again, inserting ∂C2 ∂τ and ∂ 2 C2 ∂ε 2 computed in Eq. (S.7) and Eq.(S.8), respectively, into Eq.(S.6), allows to rewrite the reaction-diffusion equation in C 3 as with the following conditions: • the initial condition: • the source: • the Neumann boundary conditions: It is preferable to have the initial condition different from 0 and the boundary conditions equal to 0. Hence, we define g(ε) and f (τ ) as: Then, Eq. (S.9) becomes with the following conditions • the initial condition: • the source: • the Neumann boundary conditions: The solution to systems of the type of Eq. (S.10) can be found in 3 .The solution to the problem will be Following what the authors have done in ceccarelli2021 we find that Hence Then, since Finally, using (S.2) and the change of variables x = L − z, we write the concentration C 1 as (S.11)

Numerical validation
In order to solve numerically the reaction-diffusion equation ( Eq. (6) of the STAR Methods section of the main text), we discretize the continuous signal density ρ S (x, t) using the Smoothed Particle Hydrodynamic (SPH) method.As explained in the STAR Methods section of the main text, we passed to an agent-based model by approximating the continuous density ρ S (x, t) with a set of discrete agents called signaling particles (SPs) that are continuously produced at the tip of the tissue, diffuse through the tissue with speed D and are removed randomly from the domain according to a Poisson process of rate k.
Considering a finite domain consisting of non-proliferating ependymal cells, we compare the 1D density profile of SPs with the analytical solution given by Eq. (S.11).In Figure S1 D, E, F, we observe a good agreement at different times between the numerical and analytical solution of the reaction-diffusion model in finite domain.This allows to validate the numerical scheme for the reaction-diffusion process of the SPs.

Figure S1 :
Figure S1: Validation of the spatial scale of the model, Related to Figure 2 and Section 2.1 of the main text (A,B) Scheme showcasing two cases where ependymal cell i (cyan disc) is a front cell (panel A) or not (panel B), using the identification rule described in section S2.1.(C) Comparison between the outgrowth predicted by our simplified 1D model (red shadowed area) and by 1 (cyan error bar).Each prediction has been obtained averaging over 100 replicates.Means are represented as lines, shadowed area and error bar correspond to 1 standard deviation.(D,E,F) Comparison between the 1D density profile of signalling particles (both the mean -blue solid line -and the standard deviation -blue shaded area -are calculated from the y-axis at each x position) with the analytical solution C 1 (z, t) (Eq.(S.11)) of the reaction-diffusion model in finite domain (orange solid line).For D = 0.6 µm 2 • s −1 , k = 1 days −1 , at different times (from top to bottom: day 0.02, day 0.2, day 4).

Figure S2 :
Figure S2: Role of the cell sensitivity to the signal density, Related to Figure 3 and Section 2.3 of the main text.(A) Cell activation function as function of the local sensitivity of the ependymal cells, for different sensitivity levels to the signal (parameter ρ min s).The maximal signal density (imposed on the right boundary, i.e. on the spinal cord front) is indicated in red.For small sensitivity to the signal (large ρ min s , yellow curve), no cells will be recruited by the signal, while for large sensitivity (small ρ min s , blue curve), the majority of cells exposed to the signal will be recruited instantaneously.(B) Error between the model predicted recruitment limit and experimental recruitment limit curve [5] as function of the sensitivity to the signal ρ min s (x-axis) and the diffusion coefficient D (y-axis).

Figure S3 :
Figure S3: Theoretical prediction of tissue outgrowth, Related to Figure 4 and Section 2.4 of the main text.Tissue domain occupied by recruited fast proliferating cells (L F (0)), tissue domain occupied by non-recruited slow proliferating cells (L S (0)) and recruitment limit (ξ(0)) at time 0. The x-axis represents the AP axis, the blue curve is the density of the signal.The amputation plane is indicated by a dashed black line.The signal sensing threshold (ρ min s ) for cell recruitment is indicated in cyan.

Figure S4 :
Figure S4: Signal-dependent transient modulation of cell cycling, Related to Figure 5 and Section 2.5 of the main text.(A, C) Percentage of cells along the AP-axis that are in G0/G1 phase (green curves) and in S/G2 phase (magenta) as function of the AP position for different times (different subplots).Dotted lines are experimental data from 1 , solid lines are simulations for different values of D and fixed ν d = 0.1days −1 (A) or different values of ν d and fixed D = 0.08µm 2 .s−1 (C).For each set of parameters, simulated data are averaged over 10 realisations (solid lines represent the mean, errorbars the standard deviation).(B, D) Relative error (see section 8 for details on its computation) between the simulated and experimental cell distributions as function of time post-amputation for different values of D while fixing ν d = 0.1days −1 (B), and for different values of ν d while fixing D = 0.08µm 2 .s−1 (D).(B', D') Number of newly recruited cells by the signal computed every 0.2 days as a function of time for different values of D while fixing ν d = 0.1days −1 (B'), and for different values of ν d while fixing D = 0.08µm 2 .s−1 (D').

Table S1 :
Model parameters, Related to STAR method section of the main text (section Mathematical modelling).Model parameters and experimental values.