Neuropathic pain caused by miswiring and abnormal end organ targeting

Nerve injury leads to chronic pain and exaggerated sensitivity to gentle touch (allodynia) as well as a loss of sensation in the areas in which injured and non-injured nerves come together1–3. The mechanisms that disambiguate these mixed and paradoxical symptoms are unknown. Here we longitudinally and non-invasively imaged genetically labelled populations of fibres that sense noxious stimuli (nociceptors) and gentle touch (low-threshold afferents) peripherally in the skin for longer than 10 months after nerve injury, while simultaneously tracking pain-related behaviour in the same mice. Fully denervated areas of skin initially lost sensation, gradually recovered normal sensitivity and developed marked allodynia and aversion to gentle touch several months after injury. This reinnervation-induced neuropathic pain involved nociceptors that sprouted into denervated territories precisely reproducing the initial pattern of innervation, were guided by blood vessels and showed irregular terminal connectivity in the skin and lowered activation thresholds mimicking low-threshold afferents. By contrast, low-threshold afferents—which normally mediate touch sensation as well as allodynia in intact nerve territories after injury4–7—did not reinnervate, leading to an aberrant innervation of tactile end organs such as Meissner corpuscles with nociceptors alone. Genetic ablation of nociceptors fully abrogated reinnervation allodynia. Our results thus reveal the emergence of a form of chronic neuropathic pain that is driven by structural plasticity, abnormal terminal connectivity and malfunction of nociceptors during reinnervation, and provide a mechanistic framework for the paradoxical sensory manifestations that are observed clinically and can impose a heavy burden on patients.


Supplementary video 4: Reinnervation of denervated territory with nociceptors.
Recovery GFP-labelled nociceptors at baseline and at 42 weeks after sham surgery over the tibial territory in SNS-mGFP mice, showing lack of overt changes in fiber patterning. Although there have been reports over the recent years of Nav1.8 expression in a few other cell types in the periphery, such as keratinocytes, this does not hold true for SNS-Cre mice, likely because unlike genetically engineered knock-in lines, the SNS-Cre line is a transgene and therefore does not include all of the promoter and enhancer elements of the endogenous Scn10a gene. Here, in SNS-tdTomato mice, we observed that 46 % of tdTomato-labelled neurons co-expressed Calcitonin gene regulated peptide (CGRP) or Substance P (SP) (putative peptidergic nociceptive cells), and another 50.5% of tdTomato-labelled DRG neurons demonstrated IB4binding (putative, non-peptidergic nociceptive neurons), consistent with the recombination profile that we have described previously (Agarwal et al. 2004;Gangadharan et al, 2011;Luo et al., 2012). These populations entail both Cnociceptors as well as A-delta nociceptors. A subpopulation of Cre-expressing neurons also express tyrosine hydroxylase (TH), which is known to be expressed in sympathetic neurons and in a small number of non-peptidergic small and mediumdiameter neurons of the DRG, which are believed to be largely C-low threshold mechanoceptors (C-LTMRs) that can be activated by light touch. In the murine lumbar DRGs, the percentage of TH-positive neurons has been reported to range from 5-15% across different studies. Here, we found that about 8% of dtTomatopositive neurons were positive for TH, consistent with the C-LTMR population being a subpopulation of C-fiber neurons (Delfini et al., 2013;Usoskin et al., 2015).
Even in pathological states, such as after nerve injury, the reporter expression in SNS-Cre mice is retained and does not appear ectopically in other sensory types.
Nevertheless, to control for this potential scenario, we also verified the key results derived from mGFP imaging experiments in SNS-mGFP mice with immunohistochemical staining of endogenous marker proteins that are expressed in specific sensory populations and employed multiple markers to consolidate the data (Fig. 3 and Extended Data Figs. 7,8 and 9).
Thy1-YFP reporter mice: There are several transgenic lines expressing YFP under a Thy1 promoter fragment, which typically imparts mosaic expression in several neural structures and outside of the nervous system. Here, we employed a line which expresses YFP in 13% of DRG neurons. At least 95% of YFP-expressing neurons are NF200-positive, indicating that they are myelinated, large diameter neurons. Although 17% of YFP neurons showed co-immunofluorescence for CGRP, there was 0% overlap with SP, another marker of peptidergic nociceptors, 0% overlap with IB4-binding, indicating absence of YFP expression in a majority of peptidergic nociceptors and all of non-peptidergic nociceptors. Consistent with a minor overlap with tdTomato expression in SNS-tdTomato mice (2%), a few thin YFP-positive fibers were seen in deeper parts of the paw skin, which also showed recovery in form of collateral sprouting from 30 weeks onwards post-SNI in the denervated tibial territory (Extended Data Fig. 3). THexpressing DRG neurons (presumably C-LTMRs) also completely lacked YFP expression. These findings at the level of DRG somata, taken together with the highly distinct morphology of large-diameter YFP-expressing fibers and their terminations at the dermis-epidermis border in the paw skin, indicate that an overwhelming majority of the YFP-positive fibers are Aß-LTMRs. YFP-positive fibers showed the classical pattern of innervation of Meissner corpuscles in the glabrous skin, in which they were found to wind around the S100-positive glial cells that comprise the corpuscle (upper panels in Extended Data Fig. 7c).
Because reporter expression can change in transgenic lines in models of injury, we backed up the data obtained with YFP imaging with immunohistochemical staining of NF200 as an endogenous marker of Aß-LTMRs in mice post-injury (Fig. 3d,Extended Data Fig. 4a) and also performed electrophysiological recordings to identify Aß-LTMRs functionally in sham and SNI mice (Fig. 3f,g;Extended Data Fig. 10).
Interestingly, we also observed that when YFP-positive Aß-LTMRs degenerate and disappear in the tibial innervation territory within 3 days of ligation and cutting of the tibial nerve in the SNI model, new YFP expressing-structures appeared in the denervated territory. These did not comprise any blebs of nerve that persist following Wallerian degeneration, but rather resembled clear cellular structures with a stellate morphology (Extended Data Fig. 3). A literature search revealed that this morphology is typical to dendritic cells, an immune cell scavenger population that are associated with clearing debris after cell injury (Hsieh and Lin, 1999;Stankovic et al., 1999;Doss and Smith, 2014). At present, the significance of these cells to nerve regeneration and sprouting or nociceptive sensitivity in the context of the present study is unclear, but will be interesting to analyse in future studies.

I. Automated stitching of four 16-bit raw image stacks into a superstack
As illustrated in Supplementary Fig. 1, a maximum intensity projection (MIP) image is calculated for each stack and SIFT/SURF is applied for pairs of MIPs. Using MIPs was a precondition to automate the procedure, because the border area would always include a sufficient number of structures for finding correspondences. The parameters identified in pair-wise MIP stitching were then applied to the 3D stacks to obtain a 16bit "superstack" covering the entire end-phalanx. To avoid disturbances by the original image overlap boundary, we applied a stitching algorithm using watershed and graph cut (Gracias, 2009, Boykov et al., 2001, Kolmogorov, 2004 to reduce the noisy influences of original image boundary. In this stitching and blending workflow, pairwise image blending is performed independently in each region by means of watershed segmentation and graph cut optimization.

II. Segmentation and removal of epidermal autofluorescence via weaklysupervised learning and label propagation
Our workflow removes epidermal autofluorescence and hair, yet retains nerve fibers in each 16-bit superstack. A weakly-supervised regularization approach was developed on discrete graph spaces for perceptual image segmentation. The weakly-supervised learning combines the labeled information (mask of surrounding epidermal autofluorescence and nerve fibers) and unlabeled data through label propagation for perceptual image segmentation. In this step, each image in a superstack will be segmented according to the segmentation parameters and label information obtained in 3D space. Epidermal irregularities and hair in the 3D image stack can be accurately classified and removed (Supplementary Fig. 2A-C).
Top-down and bottom-up segmentation strategy: statistical label propagation To improve the segmentation results, we added global 3D label information for enhancing the label propagation (Fan et al., 2018, Levin et al., 2008, Zhu et al., 2002 with global constraints in a semantic way. In this approach, a spectral clustering method is embedded and extended into regularization on discrete graph spaces. Thus, te spectral graph clustering is optimized and smoothened by integrating top-down and bottom-up processes via semi-supervised learning. Second, a nonlinear diffusion filter is used to maintain semi-supervised learning, labeling and differences between foreground or background regions. Furthermore, the spectral segmentation is penalized and adjusted using labeling prior and optimal window-based affinity functions in a regularization framework on discrete graph spaces.

Notations and their corresponding descriptions
An image is formulated as a weighted undirected graph = ( , , ) with nodes ∈ , and edges ∈ ⊆ × . Each node ! represents an image pixel ! . An edge !" connects two nodes ! and " in 8-connected neighborhood system. A is the associated affinity matrix.
is symmetric and non-negative. The diagonal matrix = ( ## … , $$ ) is called the degree matrix of graph , where !! = ∑ !" % "&# , Then = − is called the un-normalized graph Laplacian of . Assuming is connected (i.e., any node is reachable from any other node), L has certain properties of a connected graph such as symmetric and positive semi-definite, etc.
To remove epidermal autofluorescence and hair, one input image from a input 3D stack including a foreground image and a background image . The intensity of the th pixel is assumed to be a linear combination of the corresponding foreground and background intensity value, where ! is the pixel's foreground opacity. In this equation, all quantities on the right side of the compositing equation are unknown and some assumptions on the nature of , and/or are needed.

Notations and their corresponding descriptions
An image is formulated as a weighted undirected graph = ( , , ) with nodes ∈ , and edges ∈ ⊆ × . Each node ! represents an image pixel ! . An edge !" connects two nodes ! and " in 8-connected neighborhood system. A is the associated affinity matrix. is symmetric and non-negative. The diagonal matrix =

, Then
= − is called the un-normalized graph Laplacian of . Assuming is connected (i.e., any node is reachable from any other node), L has certain properties of a connected graph such as symmetric and positive semi-definite, etc.
To remove epidermal autofluorescence and hair, one input image from a input 3D stack including a foreground image and a background image . The intensity of the th pixel is assumed to be a linear combination of the corresponding foreground and background intensity value, where ! is the pixel's foreground opacity. In this equation, all quantities on the right side of the compositing equation are unknown and some assumptions on the nature of , and/or are needed.

Regularized semi-supervised cost function
To derive our solution for the grayscale case we make the assumption that both and are approximately constant over a small window around each pixel. Note that assuming and are locally smooth does not mean that the input image is locally smooth, since discontinuities in can account for the discontinuities in . This assumption expresses as a linear function of the image , ! ≈ ! +b, ∀ ∈ , '() and is a small window. The linear relation is similar to the prior used in (Grady et al., Levin, 2008, Fan et al., 2018. Our goal is to find and minimizing the cost function where " is a small window around pixel . This cost function has a regularized term on for numerical stability. This function can be converted as a quadratic cost function in the formula with unknowns, ( ) = . . Here we can consider a Laplacian diffusion matrix , which is a symmetrically positive definite matrix.

User label constraints
To achieve reasonable segmentation of epidermal autofluorescence and appendages and internal nerve fibers, biologists mark and label several constraints in image 3D stacks. We mark autofluorescence components in black pixels ( = 0) and consider neural fibers as foreground ( = 1). Based on the utilization of label constraints, we solve the following equation, where is chosen as a number for optimizing the segmentation region, / is a diagonal matrix whose diagonal elements are one for constrained pixels and zero for all other pixels, and / is the vector containing the specified alpha values for the constrained pixels and zero for all other pixels. Since the above cost is quadratic in alpha, the global minimum may be found by differentiation of the equation and setting the derivatives to zero. This amounts to solving the following sparse linear system, In this equation, let be an image formed from and according to the compositing equation (1), and let * denote the true segmentation. If and satisfy the model in every local window 1 and if the label constraints are consistent with * , then * is an optimal segmentation solution for the system (3), where is constructed with = 0.
For example, the extracted segmentation region and the smallest eigenvectors tend to be piecewise constant over the same regions. If the values inside a segment in the eigenvector image are coherent, a single label within such a segment should suffice to propagate the desired value to the entire segment. On the other hand, areas where the eigenvector's values are less coherent correspond to more "difficult" regions in the image, suggesting that more label efforts might be required there.

Spectral Laplacian affinity clustering
The matrix = − is the graph Laplacian for segmentation. The matrix has the same form as the graph Laplacian used in spectral methods for segmentation but with a novel affinity function given by (Fan et al., 2018). For the typical way to define the affinity function (for example, for image segmentation using normalized cuts (Shi et al, where is a global constant. This affinity is large for nearby pixels with similar intensity and approaches zero when the gray value difference is much greater than .
The transformation is image-dependent and is estimated using a manifold learning technique. Based on the Laplacian as = − , we obtain the affinity function with certain constraints provided by label propagated information.
These two affinity functions 2 and 6 can be examined by using eigenvectors as guides for achieving optimal segmentation. The smallest eigenvectors of the Laplacian can guide the user where to label and mark.

Iterative optimization with soft constraints in multi-resolution
The optimization problem defined by equation (3) minimizes a quadratic regularization function subject to linear constraints, and the solution can therefore be found by solving a sparse set of linear equations. Large images or image 3D stacks cannot be solved using Matlab's solver due to memory limitations. To overcome this, we use a coarseto-fine multi-resolution scheme. We down-sample the image and the constraints and solve it at a lower resolution. The recovered region is then interpolated to the finer resolution. The values are thresholded and pixels with close to zero or one are clamped and considered as constrained in the finer resolution. Therefore, we minimize the regularization function in the equation where != , !> , != , !> , are the and derivatives of and , and != and !> are the segmentation region derivatives. We note that for a fixed , this function is quadratic, and its minimum can be found by solving a sparse set of linear equations. Given the solution of and , the solution be further refined according to the practical problems.
For certain complex distribution of phalangeal nerve fibers, epidermal autofluorescence and appendages, such a segmentation-based reconstruction may produce noisy results, and therefore, we solve for and using the compositing equation with some explicit soft and hard constraints. Soft constraints include smoothness priors on and . The reconstructed image is composed of foreground neural fibers in the central region and background distributed in the surrounding regions. Constrained pixels maybe eliminated from the system, reducing the system size. For that, we note that within the constrained areas, there is no need to enforce the local linear models.
Therefore, when computing the segmentation Laplacian matrix in the equation (4) and (6), we sum only windows 1 that contain at least one unconstrained pixel. Aside from efficiency, clamping alpha values to zero or one is also useful in avoiding an oversmoothed alpha region. We note, however, that clamping the alpha values in a coarse to-fine approach has the adverse effect that long thin structures such as strains of hair may be lost.
Improved 3D segmentation using Z-direction geodesic coherent label propagation in image 3D stacks Based on our datasets, we needed to achieve maximal accuracy to segment and remove the epidermis with strong intensity, partial epidermis inner buffer region, hair and other autofluorescence. Therefore, stronger prior information (e.g., hard constraints) were combined with soft constraints such as smooth priors for computing an accurate segment boundary between foreground and background.
Our datasets are multiple image 3D stacks. In practice, we give certain sparse labels from the top of the scanning layers and also some labels at the bottom layers in 3D image stacks. The goal of this label information is not only for 2D segmentation between foreground and background in 2D images, but also for utilizing z-direction geodesic coherent information ( Supplementary Fig. 2D-F), pairwise similarity and connectivity features. In z-direction following the depth of image 3D stacks, the geodesic metric data similarity is much more effective than previous Euclidean distance-based similarities.
The computation of eigenvectors (Chennubhotla, 2005) of the Laplacian depend only on the input image and are independent of the user's label constraints. It is not necessary to compute them, unless the user wishes to use them for guidance in scribble placement, as described earlier. In this case, they only need to be computed once, possibly as part of the initialization that takes place when a new image is loaded. With geodesic zdirection propagation in image 3D stacks, the coherent pairwise similarities propagate a closing boundary in z-direction from the top layer to the bottom layer including epidermal autofluorescence and appendages.

III. 16-bit to 8-bit conversion of superstacks by nonlinear adaptive depthdependent adjustment on global and local scales
To facilitate the analysis steps described in sections 4 and 5, the data was converted from 16-bit to 8-bit TIFF files using an approach that maximally preserved the dynamic range of the signal. We proposed a global and local contrast adjustment and combinatorial optimization approach for adaptively adjusting the non-uniform illumination image stacks. The balance of local contrast adjustment and the global contrast adjustment is optimized based on a contrast-brightness based pixel-level weighted fusion framework without losing any relevant information. The optimization framework is formulated in the following, Where -( ) is an optimal enhanced image, is the index of pixels in a given 16-bit image. We compute the weight maps for the intensity images in a global weight _ 2 based on a global image intensity -( ), and the local weight map _ ? based on a local image intensity -( ) . The result in this equation is an optimal fused image in a Laplacian pyramid method (Burt, et al., 1983).
Step 1: Global adaptive contrast and brightness enhancement for gray images A general histogram modification framework (Arici, 2009) was used to obtain a mapping function for the optimal modified histogram. We consider a gray image , with a total number of pixels and an intensity level range of [0, − 1], R is a 16-bit value 65536. The normalized histogram ℎof the image is computed in equation (2).
where ∈ [0, − 1] is pixel values, and @ is the total number of pixels having the same value s. The modified histogram ℎ should be closer to the normalized uniform histogram ℎ A and the value of the residual ℎ − ℎshould be small. The problem of obtaining the optimal modified histogram ℎ h can be regarded as a solution of a bi-criteria optimization problem, which can be formulated as a weighted sum of two terms as Eq. (3).
where ‖ℎ − ℎ -‖ is a norm of ℎ − ℎ -, ‖ℎ − ℎ A ‖ is a norm of ℎ − ℎ A , and the parameter λ > 0 adjusts the trade-off between the contrast enhancement and the fidelity of the data.
Step 2: Local contrast adaptive adjustment and global optimization framework Local contrast was enhanced using the framework represented in equation (4) and the Contrast-Limited Adaptive Histogram Equalization method (CLAHE) (Pizer, 1987).

IV. Automated rigid and non-rigid 4D registration of stacks acquired at different time points
First, automated registration was implemented in a stratified 3D model registration framework, which efficiently handles a hierarchical pyramid multi-resolution 3D image stack. Our algorithms compute the multimodal 3D registration in a stratified approach with respect to the rigid and non-rigid 3D alignment including 3D-to-3D pose correction, 3D-to-3D rigid transformation and 3D-to-3D linear and nonlinear registration (Supplementary Fig. 3A-D).
• First, the 3D pose of a given image stack is automatically corrected using the relative 3D pose of the consecutive image stack.
• Second, the rigid linear alignment in the iterative closest points (ICP) (Bergevin, 1996, Besl, 1992 and the nonlinear registration in the thin-plate-splines (TPS) (Bookstein, 1989, Szeliski, 1996, Zheng et al., 2009) are stratified into two steps for reducing and avoiding the over-blending effects in linear and nonlinear 3D registration.
• The results of the stratified matching engines are integrated to achieve higher accuracy of 3D registration of nerve fibers.
Following this strategy, the original global distribution and locally detailed structure of nerve fibers are largely and optimally preserved after 3D registration. The approach performs efficiently in particular for large-scale 3D data alignment.
Rigid 3D registration for 3D image stacks acquired at different time points To find the largest similarity (e.g., intensity localization and patterns) between two 3D stacks, representative feature points were extracted from the source and target 3D image stacks in 3D space. To achieve this, we firstly extracted 3D point clouds of nerve fibers from converted 8-bit high-resolution 3D image stacks. We then estimated the 3D registration parameters from the extracted 3D point clouds based on the iterative closest point algorithm described below. According to the ICP algorithm, the rigid pairwise local surface registration algorithm can be described as follows. The 3D registration method seeks to find the best transformation T that relates two entities P and Q whose 3D point clouds are given by F and G , respectively. The goal is to find Τ F such that the objective function D F , G E is minimized, where Ψ: P → Q; for ⋁ ∈ F , Ψ(P) ∈ G . The transformation Τ F is used to optimally align two point clouds P → Q. In practice, the function Ψ(P) is unknown and . While the mean square error is smaller than a given threshold, the iterative optimization converges monotonically to a global minimum. When a good initial value is given, the algorithm can get a global convergence for aligning pairwise 3D point clouds or image 3D stacks.

Non-rigid 3D registration based on Rigid ICP registration
We use the thin-plate-spline (TPS) method for global registration by mapping each local feature point onto its global position. The TPS method has been introduced to 3D geometric deformation processing (Bookstein, 1989, Pottmann, 2006. The TPS has an algebra expressing the dependence of the physical bending energy of a thin plate on The algorithm of TPS is a non-rigid, globally smooth function including affine and nonaffine components. The non-affine warping component means that the sum of squares of all second order partial derivatives is minimized. Such functionality can be used for 3D mapping and it is easily computable. Thus, if : $ → is an n-dimensional TPS, the bending energy is minimized described in the following,  (Duchon, 1977), This equation means that the TPS fits a mapping function ( ! ) between corresponding point-sets and by minimizing the energy. The first term is a distance measure term between two sets of points. Since ( ! ) is not unique, we need smoothness constraints to regulate the mapping. Here, the spline will not be interpolated, but for any fixed , there is a unique minimum. The warping strength depends on the regularization parameter . While is close to zero, we get an exact alignment of corresponding surface vertices. If is zero, interpolation is exact and as it approaches infinity, the resulting TPS surface is reduced to a least squares fitted plane, i.e.,"bending energy"of a plane is 0. For the interpolating case, the thin-plate spline specification provides a linear system of equations.
To apply this equation to solve 3D point sets, based on (Wahba, 1990), we performed a QR factorization to decompose the matrix on Ψ using the approximated equation, with the condition of . = 0. Where is × and !" = (| ! − " |) is a Green's function. We set a small value of to improve the stability of numerical computation for this function.
In contrast to previous work, the alignment of numerous local surfaces is globally controlled by the low-resolution global shape model. Since we need to have precise alignments for rigid local scans using the feature prior from the low-resolution global shape model, the spline must be heavily weighted towards interpolation for connecting and aligning slightly-overlapping high-resolution local surface patches. For this reason, we can rely on accurate correspondences to produce a good local surface alignment without global shape distortion or error.
Estimating transformation parameters based on extracted representative feature point clouds In practice, since one 16-bit image stack consumes 16 GB, an entire end phalanx superstack is around 60 GB. Hence, when aligning superstacks acquired at two timepoints, a computer RAM size exceeding 128 GB is required. Furthermore, computing the transformation parameters directly from the superstacks is time-consuming.
Therefore, we first perform SIFT or SURF feature detectors in 3D spaces to get representative feature points with only 3D coordinates in 3D spaces. The ICP algorithm can be performed directly using the extracted feature coordinates, e.g. one 3D point cloud is the target point cloud and the other one is a movable point cloud ( Supplementary Fig. 3A-D). In this way, 3D transformation parameters can be determined in an efficient way and applied to pairwise superstacks only for pixel-based transformation.
Rigid global 3D alignment and nonrigid local elastic adjustment for optimal 3D registration In the work diagram, firstly, high-resolution local surfaces are aligned onto the lowresolution global shape model using the rigid ICP method. In this way, the original fidelity of high-resolution local surfaces can be kept after the alignment. Secondly, due to the memory bottleneck problem, we extract representative feature point clouds to reduce the computation load and apply the ICP transformation parameters directly to original image 3D stacks. After the ICP alignment for all pixel-based 3D transformations in the first step, a good starting position for non-rigid TPS alignment in the second step is initialized. In this way, the non-rigid alignment in TPS will be optimized for local adjustment of alignment, shown in Supplementary Fig. 3E-Q.
Simultaneously, each aligned local nerve fiber and fiber endings are non-rigidly adjusted and refined using a global non-rigid TPS bending algorithm based on initial ICP alignment. 3D alignment from rigid registration to non-rigid 3D registration, for morphology deformation measure and quantitative analysis in 3D and 4D space.
Robustness of 3D registration algorithm 2-photon end phalanx in vivo imaging implies changes between consecutive imaging sessions. In practice, due to the large space of possible 3D nerve fibers changes and deformations, 3D registration of nerve fibers with such partial similarity tends to be challenging, as illustrated in Supplementary Fig. 4. Therefore, we had to find and utilize more representative local features and global topological-structure information for 3D registration. Based on the RANSAC algorithm (Fischler et al., 1981), we further optimized the ICP rigid alignment with RANSAC matching and correspondence searching criteria. The alignment can then be obtained by solving a global convex optimization problem with local geometric deformation measure using interactive label information.
In addition, to improve the correspondence searching between two 3D stack pairs, the nerve fibers are decomposed into adaptive patches according to the tree-like structure of the network formed by the nerve fibers. Secondly, the approach is capable of simultaneously computing the mean distribution of fibers in searching 3D space, represented by a probability density function, from multiple fiber point sets (represented by finite-mixture models), and registering them non-rigidly to this emerging mean distribution shape. Third, the normalized 3D pose and correspondence is served as an initial correspondence for 3D image stack pairs. Since we want to know where and how much of the changes happen, we need to align these 3D datasets using existing partial similarity features including global topology information and local feature similarity information with Z-depth information (Supplementary Fig. 4). The approach is also capable of handling large changes of nerve fibers, in which each test data is represented as compact training features. In addition, it is faster than stateof-the-art approaches and achieves high-quality registrations of two neural fibers with growing changes and with different acquisition angles.
Estimation of maximum common volume (MCV) for image 3D stacks acquired at multiple time-steps in vivo Shifts in the imaging volume can occur with every imaging session, posing the challenge to crop all superstacks to the same common overlapping volume in order to prevent biases in subsequent quantitative analyses. To solve this problem, we combined two registration workflows to find the best and maximum common regions in each superstack of a time series. We combined the all-to-one 3D registration and pairwise 3D registration to find the MCV (Supplementary Fig. 5A). Subsequently, changes in nerve fibers were certain to arise from the same MCV, so changes could not arise artificially from non-overlapping volumes.

V. Automated neuronal segmentation, tracing and statistical neural network analysis
Automated and accurate neural fiber tracing and analysis Many regularization methods work by pruning noisy branches of the medial axis (Li et al. 2015;Yan et al. 2016), which requires an initial approximation of the medial axis.
Based on our datasets and previous studies (Amenta, 2001, Andrea et al, 2016, Huang et al, 2016, Siddiqi and Pizer 2008Tagliasacchi et al. 2016), we developed a datadriven and task-oriented nerve fiber 3D tracing approach.
Our method rests on two insights on neuronal fiber tracing. Firstly, we show that the medial axis representing shape and structure of a neural fiber can be well approximated, both topologically and geometrically. Second, we show that our approach can provide not only a topologically correct and geometrically convergent approximation from global-to-local optimization, but also magnifies the detection sensitivity of weak nerve fiber signals in deep scanning layers. It is also one of the main differences between our approach and state-of-the-arts approaches in that we can largely trace and analyze discontinuous or broken fibers in the end phalanx (Supplementary Fig. 5B-G). Hausdorff distance. To measure the distance between two sets # , * , we use the symmetric Hausdorff distance d S defined as, Chazal (Chazal et al., 2005) show that, even though the medial axis is highly sensitive to boundary perturbations, a subset of the medial axis enjoys certain stability properties when the perturbation is bounded by the Hausdorff distance.
Incremental and ensemble L # -regualrization for accurate 3D nerve fiber tracing We extended L # -regularization for computing 3D skeletons of nerve fibers and also magnifying the tracing capability for weak fibers in the deep layers. Given an unorganized and unoriented set of points V = { W } ⊂ < , ∈ , we construct the following definition for incremental and ensemble computing the optimal 3D medial skeletons that leads to an optimal set of projection points = { X }, ∈ , Enhancement of medial axis detection in 3D spaces using ensemble tracing prior Where the first term is a localized median points of V, the second term ( ) regularizes the local multiple point distribution of , indexes the set of projected points , and indexes the set of input points V. The weight function ( ) = (J # (Z * ⁄ ) # ⁄ is a smooth function with support radius ℎ defining the size of the supporting local neighborhood for incrementally increasing the single tracing median skeleton in a L # -regularization iterative optimization way.
The main advantage of definition (3) is that it does not place strong requirements on the quality of the set of points representing the input fiber points. In particular, it can be directly applied to raw multiple tracing point clouds with a scanning depth of 700 layers (e.g., 700 µm) and our approach leads to magnifying the detection sensitivity of neural fibers in a 3D skeleton.
Applying the equation (3) tends to yield a sparse distribution, where the local centers are accumulated into a set of point clusters. To avoid such clustering and maintain the proper medial geometry representation, further accumulation needs to be prevented once points are already contracted onto their local center positions. This is achieved using ( ) in (1), which adds a repulsion force whenever a skeleton branch is formed locally.
For detecting and refining the branches of neural fibers in skeleton computing, we adopted classical weighted PCA to detect the formation of skeleton branches and extended the method of mean shift (Fukunaga andHostetler 1975, Grillenzoni, 2018) into 3D space. Mean shift is an iterative procedure that shifts each data point to the weighted average of data points in its neighborhood.
Iterative reconstruction of the broken skeleton in broken fibers with weak signals is specially developed for our datasets. We first use neighborhood data points to estimate a representative central median skeleton in a weighted mean shift adjustment in our regularization function. We also apply a simple 3D buffer to guide the ensemble extension of the median skeleton in broken region or volume as illustrated in Supplementary Fig. 6. Inside the blue buffer volume in Supplementary Fig. 6 H-L, the ensemble growing and mean-shift will connect the broken skeleton by using weak fiber signals. These weak fiber signals will give strong weight inside the buffer volume for connecting the broken in a locally-weighted ensemble construction way. Following the buffer volume, the computation will incrementally connect the broken fibers with weak fiber signals in a shift window to the global volume. As a result, our algorithm may produce skeletons with certain global and local constraints with high-fidelity and correct topology as highlighted in Supplementary   Fig. 5E-G. We identify the number of traced fiber clusters during the growth period.
The geodesic distance measure can largely help to measure the broken distance in a broken volume which is often by the constraints of utilizing the buffer volume. In this way, we can also add more information for optimal reconstruction of the broken fibers, e.g., normal of fiber surface patches in Supplementary Fig. 6, and outward directions directly derived from spatial-temporal fiber changes.

Construction of 3D graph for spatially-localized and topological network analysis
To achieve cluster analyses of 3D fiber networks, localization of fiber changes, as well as activity-dependent structural plasticity analysis were developed. These algorithms include: • Analysis of connectivity and fibers using graph G=(V,E), e.g., dendrogram, and make the graph more bio-informative, • Accurate 3D registration in a sliding 3D window to measure fine structures and branching changes in 4D space.
• Localize and measure dynamic changes of fine structure of dendritic filopodia and small branches in 4D space.
• Dynamic change of neural pattern and neural connections (clusters, undirected graph, data-driven and task-oriented quantitative measure).
A graph network-based cluster analysis in 3D space revealed the location and quantity of broken fiber networks, location and quantity of regenerated fiber networks, etc. Furthermore, we construct the undirected graph network to represent the neural fibers in a graph network. In this structural graph network, we can measure and localize the changes not only using Euclidean distance measure (e.g., 3D registration using Euclidean distance) but also using geodesic distance measure of fiber length, fiber connectivity and localize the fine structural changes in a converted graph network.

Construction of 3D graph network
In this step, the traced 3D skeleton was converted into a graph network to measure neural fiber plasticity and changes in 4D space. The network topology is represented as a sparse connectivity matrix , where !," is the length of the shortest link between Global connectivity and local connectivity measured in graph network using geodesic measure.
Here we use the geodesic path as the shortest path and distance along the fibers. The shortest path (SP) and the average shortest path (ASP), which measure the typical separation between two nodes in the network, is derived by taking the mean of all shortest paths between any two nodes. The topological architecture of the network allows quantifying the efficiency of the network for information transfer based mainly on geodesic measures. We found that the shortest paths between nodes in the network run along the fibers with most high-centrality nodes clustered together. The centrality of a connected node is the fraction of all such shortest paths in the network running through that node. This analysis provides insights into neural plasticity changes using local to global principles into a large, interconnected network during neural degeneration and regeneration. VI. Quantitative analysis of structural plasticity in SNS-mGFP and Thy1-YFP mice Differences of nerve fiber morphology between SNS-mGFP and Thy1-YFP mouse lines are illustrated in Fig. 1. In the 3D spatial domain, these differences mainly include the distribution density of fibers, length, thickness of fibers and the location of fibers.
Therefore, the segmentation and tracing algorithms were adjusted accordingly as detailed in the following.

4D analysis of epidermal fibers in SNS-mGFP datasets
The analysis workflow designed to capture changes of thin epidermal fibers is described in the following (Supplementary Fig. 7). Epidermal fibers and deep fibers in SNS-mGFP mice paw are classified so that changes of epidermal fibers can be quantitatively measured and localized. The workflow applied to a complete SNS-mGFP dataset acquired at all different timesteps.

Supplementary note 3: Tyrosine hydroxylase (TH)-expressing fibers
TH is expressed in all sympathetic and dopaminergic nerve fibers (Brumovsky, Pain 2016). Accordingly, anti-TH immunohistochemistry in glabrous hindpaw plantar skin of sham-treated SNS-mGFP revealed their localization in deeper parts of the dermis, consistent with adrenergic innervation of sweat glands. A sub-population of mGFPpositive sensory fibers were also immunoreactive for TH in the dermis. This is consistent with the current literature which indicates that TH-expressing neurons in the DRG are C-LTMRs, which is a subset of non-peptidergic neuronal population expressing Ret and GINIP, TafA4 and Vglut3, but distinct from the IB4-binding or MGPRD population (Brumovsky, Pain 2016;Delfini et al. 2013;Piers et al. 2015;Usoskin et al. 2015). TH-positive C-LTMRs and are known to robustly innervate hair follicles in hairy skin with very little innervation of glabrous skin reported (Brumovsky, Pain 2016). Accordingly, none of the 8 sham-treated SNS-mGFP mice tested show THimmunolabelling at the dermal-epidermal border or in the epidermis of the glabrous skin (Extended Data Fig. 9a). A similar pattern of TH-expression was observed in the denervated tibial territory of 6 out of 8 SNS-mGFP mice tested at 42 weeks post-SNI (Extended Data Fig. 9c). In contrast, the remaining two mice showed presence of anti-TH immunolabelling, which was partially overlapping with mGFP expression, at the dermis-epidermis border and partly invading into the epidermis (Extended Data Fig.  9b). This raises the possibility that C-LTMRs may also collaterally sprout from hairy skin into denervated glabrous skin. Given that they are a sub-population of the SNSfibers (Nav1.8-expressing) that were ablated in the functional experiments shown in Fig. 4, it cannot be discounted that they contribute to mechanical allodynia in the tibial denervated territory. However, this is not likely given the low incidence of the putative collateral sprouting of TH fibers observed as well as the lack of functional indications thereof in the skin-nerve recordings we performed. In the absence of hair shafts, which they normally act as a lever for transmitting force to the nerve ending, it is unclear whether

Supplementary note 4: Transcriptomics
We performed RNA expression analyses to address whether nociceptors which collaterally sprout from the sural nerve into the denervated tibial territory at 24 weeks post-SNI show molecular differences from the native nociceptors of the tibial nerve which populate and innervate this territory in physiological conditions (sham). The somata of these axons were retrogradely labelled by injection of tracers in a very small volume in the plantar paw skin corresponding to the tibial territory and picked following acute dissociation. During this step, we segregated tracer-labelled neurons as being IB4positive and IB4-negative depending upon binding with fluorescently labeled-IB4 and subjected them to RNA sequencing as described under methods. As additional controls and comparisons, we also performed similar analyses on somata of fibres that innervate the sural territory in control and neuropathic mice SNI 24 weeks. All raw and annotated gene expression data with corrected p values are uploaded on a public repository European Nucleotide Archive (https://www.ebi.ac.uk/ena) under the accession number PRJEB50184. In the description below, the labels 'tibial' and 'sural' denote the corresponding skin territory. IB4-positive neurons are unequivocally identified as nociceptors in all groups. In case of IB4-negative neurons, although we tried our best to pick small-diameter cells post-dissociation, these can potentially include nonnociceptive neurons, particularly in the control condition for the tibial territory group and both control and SNI conditions in the sural territory group. However, we would like to emphasize that we have also analyzed data in the IB4-negative population and observed that there are no differences which would go against the conclusions derived from the gene expression analyses.
Our main findings are as follows: 1.
Although up to 26154 genes were studied in total, the number of regulated genes is very modest (in total 24 in tibial territory and 68 in sural territory), as shown in Extended Fig. 13. This is completely plausible, in our view, since analysis was performed at late stages (6 months post-injury), when local injury and inflammation pathology is absent and reinnervation has already occurred, so the massive regulation of injury and degradation-related genes that is typical to early stages is not expected.
Secondly, previous studies that looked at gene regulation post nerve-injury analyzed all DRG neurons in bulk. This study compared nociceptors in their native state of innervation with nociceptors that had collaterally sprouted and established reinnervation. When comparing gene expression differences, we noticed that the largest differences actually occur between IB4+ and IB4-neurons, which are given by the nature of their cell specialization, in both sham and SNI mice; they are, therefore, addressed separately in the data below. Within the same subpopulation of nociceptors, gene expression did not vary much between the treatment groups.
2. The representation of different sub-classes of nociceptors in the tibial territory post-SNI was compared with sham mice using deconvolution analyses and clustering with marker proteins described by Usoskin et al. (2015). We did not find any significant differences in the relative abundance of the NP1 and NP2 sub-populations and PEP1 and PEP2 sub-populations of peptidergic nociceptors. The PEP3 subclass was represented to a very low extent across the sham samples and not detectable in the SNI samples, but owing to the low abundance in the overall population of nociceptors, this was not significant. A Chi Square analysis of comparison of overall distribution of the 5 sub-populations of nociceptors did not reveal a significant difference (Extended Fig.   13c). This, therefore, argues against dominance of a particular class of nociceptors with altered or special properties amongst the collaterally sprouting nociceptors as a mechanistic basis for mechanical allodynia developing in the denervated sural territory.
3. Modulators of neuronal excitability: As compared to control mice, we did not observe any upregulation of genes encoding reported mediators of peripheral sensitisation in the nociceptors reinnervating the tibial territory in SNI mice (24-26 weeks post-SNI). For example, there was no upregulation of ion-channels transducing sensory stimuli, such as TRP channels, Piezo2 or ATP-gated channels (P2X3), or enhancing membrane excitability, such as diverse sodium channels, calcium channels or potassium channels that have been linked to nociceptor sensitization, amongst others. Data from IB4positive (putative non-peptidergic nociceptors) and IB4-negative populations are shown below in Table 1 and Table 2, respectively.  Moreover, in DRG neuronal somata with terminations in sural territory, we also did not find major indications of gene expression changes indicative of classical nociceptor sensitization at this chronic time point post-SNI. Genes encoding the sodium channel Nav1.1 (Scn1a) and the key growth factor receptors, TrkA (Nrtk1) and TrkB (Ntrk2), were actually significantly downregulated (rather than being upregulated) in sural territory neurons of SNI mice as compared to sham mice. We conclude therefore that nociceptor sensitization phenomena induced by molecular mediators that are hallmarks of other types of pain, such as acute and chronic inflammatory pain, are not causing the observed change in mechanical thresholds upon reinnervation.
4. In DRG neuronal somata of nociceptors that collaterally sprouted into the tibial territory at 24-26 weeks post-SNI, instead of finding regulation of genes linked to peripheral sensitization as explained above, we found that most of the genes that were significantly differentially expressed in collaterally sprouting nociceptors over sham conditions are directly linked to structural changes, such as neurite outgrowth, axonal pathfinding, neuroprotection and axonal survival, or belong to families that have been ascribed these functions. Some of these genes are indeed downstream of NGF and reported to be important for NGF to stimulate neurite outgrowth, supporting the point about structural remodelling and defects in end organ connectivity as underlying mechanisms for the observed phenotypic changes (Extended Fig. 13). There are also genes implicated in cell differentiation, cytoskeletal changes as well as cell organelles, such as lysosomes, endosomes, inflammasomes and processes that involved in cellular detoxification, such as ubiquitination and phase II metabolic enzymes. Please see complete list below for collaterally sprouted nociceptors as compared to the native innervation of the tibial territory in tables 3 and 4 below.