Physically Consistent Scar Tissue Dynamics from Scattered Set of Data: A Novel Computational Approach to Avoid the Onset of the Runge Phenomenon

: The foreign body reaction is a complex biological process leading to the insulation of implanted artiﬁcial materials through a capsule of scar tissue. In particular, in chronic implanta-tions of neural electrodes, the prediction of the scar tissue evolution is crucial to assess the implant reliability over time. Indeed, the capsule behaves like an increasing insulating barrier between electrodes and nerve ﬁbers. However, no explicit and physically based rules are available to com-putationally reproduce the capsule evolution. In addition, standard approaches to this problem (i.e., Vandermonde-based and Lagrange interpolation) fail for the onset of the Runge phenomenon. More speciﬁcally, numerical oscillations arise, thus standard procedures are only able to reproduce experimental detections while they result in non physical values for inter-interval times (i.e., times before and after experimental detections). As a consequence, in this work, a novel framework is described to model the evolution of the scar tissue thickness, avoiding the onset of the Runge phenomenon. This approach is able to provide novel approximating functions correctly reproducing experimental data ( R 2 (cid:39) 0.92) and effectively predicting inter-interval detections. In this way, the overall performances of previous approaches, based on phenomenological ﬁtting polynomials of low degree, are improved.

Several studies explored the FBR details describing a first phase related to the adsorption of blood and plasma proteins (e.g., fibrinogen, fibronectin, albumin and antibodies) to the implant surface [29,30]. Therefore, a second phase was described where the adsorbed protein layer promoted monocyte and leukocyte extravasation, attraction and adhesion to the surface, together with the coagulation cascade [31][32][33]. Then, a sort of matrix of fibrin was found to attract circulating leukocytes and local macrophages around the surface of the implanted material under the chemical attraction of chemokines [32,34]. Furthermore, it was described as the way extravasated monocytes differentiate to macrophages and fuse together to form multinucleated foreign body giant cells (FBGCs). These cells started releasing further inflammatory cytokines, boosting the inflammatory response through a positive feedback mechanism [35,36]. Then, it was described as myofibroblasts triggered a massive secretion of extracellular matrix (ECM) components [15,27,37]. Finally, at the end of the whole process, the scar tissue capsule was found to be impermeable to the nonspecific immune system and to many chemicals [14,15,27].
Some previous studies explored the evolution of the FBR caused by intraneural electrodes longitudinally implanted (e.g., LIFE electrodes [49,50]) in the rat peripheral nerves, while others explored the cellular composition and the thickness of the fibrotic capsule in a time span ranging from 1 day to 8 months [51,52]. These investigations revealed that the contacts between the nerve fibers and the electrode active sites were partially lost because of the outgrowth of the encapsulating tissue. In particular, the increase of the scar capsule was directly related to the increase of the electrical impedance, resulting first in the decrease of the signal to noise ratio [53] and possibly in a sensible drop of registration and stimulation capabilities [36,54].
Therefore, a correct and quantitative assessment of the capsule thickness over time is crucial information to improve the design of intraneural interfaces as well as their implantation technique. However, the development of the FBR involves several interrelated complex processes [14,15,[27][28][29][30][31][32][33][34][35][36][37][38] which are partially hidden since their explicit and physically based rules are not known in all details.
Nevertheless, it is a pressing issue to investigate whether computational strategies are able to reproduce the scar tissue evolution, making sense of the huge complexity of the underlying phenomena. In particular, a suitable trade off between light computational models and physically meaningful dynamics has to be found. Under these constraints the identification of a reliable procedure to reproduce the scar capsule thickness was not foreseen and standard numerical approaches [55,56] fail. In addition, if the experimental sets of data are small, with experimental time points which are closer together, ill-conditioned systems could arise, further increasing the complexity of the problem.
As a consequence, in this work, a novel numerical procedure was provided to reproduce and predict the evolution of the scar capsule mean thickness, starting from a whole set of scattered experimental data with low cardinality [57]. More specifically, two classic approaches, fully exploiting experimental data and optimizing standard metrics, are first presented for benchmark. Therefore, the onset of the Runge phenomenon, affecting their ability to reproduce data through high-degree polynomial functions [55,56], was explored. Since the optimization of standard metrics was not able to avoid huge errors in predictions due to the onset of this numerical phenomenon, a novel approach was proposed to overcome previous standard models by optimizing the parametrization of the scar thickness evolution over time. Finally, a novel metric was provided to better investigate the suitability of potential candidate functions in reproducing the experimental data set.

Methods
Experiments were performed to measure the thickness of the scar tissue capsule around polymeric stripes implanted within in vivo mice, as explained in detail in [57]. A standard procedure was used to evaluate the microstructure of the nerve and the thickness of the tissue capsule around the implant, as explained in [58]. More specifically, nerve segments were postfixed and embedded in epon resin, while semithin sections (0.5 µm thick) were stained (toluidine blue) and investigated by light microscopy. A representative section from the central part of the implant in the nerve was chosen for each animal, while the thickness of the scar capsule was computed as the distance between each side of the implant and the closest myelinated axon using ImageJ software (U.S. National Institutes of Health, Bethesda, MD, USA). In particular, the scar thickness was assessed at 1, 2 and 4 days and 4, 8, 16 and 32 weeks post insertion [57]. This experimental approach resulted in a scattered numerical subsets(t) = {s(t 0 ), . . . , s(t 8 )} with low cardinality. First, a well known standard approach (i.e., Vandermonde-based polynomial interpolation [59,60]) was used to approximate the continuous evolution of the thickness over time from scattered data. More specifically, the basis of polynomials {1, t . . . . . . , t n } was chosen and the coefficients of the interpolating polynomial function were found through the solution of the system: wheres(t) was the experimental data vector, while V(t) was the classic Vandermonde matrix defined as: while a = {a 0 , a 1 , . . . , a n } T was the vector of numeric coefficients. Therefore, once the system in Equation (1) was inverted, the interpolating function was written as: where n = 9 was the number of data points known through experiments. A second standard approach, fully exploiting the experimental information over time, was implemented through the Lagrange interpolation polynomials [59,60], which were defined as: where n was, again, the number of experimental points. Once developed, this function was written as: where x i ∈ R were numerical coefficients related through Equation (4) to s(t j ). However, since the cardinality of the known set of data points was small, supplementary points were added to improve the quality of the interpolation process. More specifically, the last two points of the experimental time detections were used to obtain a further middle time point through t ω = [t n + t n−1 ]/2 (6) where, n = 9 was the actual cardinality of the experimental data set. The thickness related to this further time point was found as: In the last part of the curve a steady state was guessed. After the first application of Equations (6) and (7) the cardinality of the set of points was increased by one unit (i.e., n = 9 + 1 = 10). Similarly, this procedure was iteratively performed two times in order to add two and three further points to the initial data set. In other words, each novel set of points was obtained from the previous one by adding one further point, so the cardinality of the novel sets was n = 9 + 2 = 11 (second iteration) and 9 + 3 = 12 (third iteration), respectively.
Since more data points were exploitable, higher degree interpolant polynomials were used. In addition, these polynomials were recursively connected through the Neville's formula: where i, j ∈ N were both i, j < n, and i = j were two different integers. In addition, g n (t) was the n-th order function, interpolating the data points at t 0 , t 1 , . . . , t n , while g n−1,i (t) was the n-1-th order approximate function, interpolating the data points at t 0 , t 1 , . . . , t i−1 , t i+1 , . . . , t n . Moreover, the previous approaches were heavily affected by the number of the experimental points which were available for the reconstruction procedure. As a consequence, a novel approach was proposed to disentangle the reconstruction of the thickness over time from the available experimental information. More specifically, novel functions JHS 0 (t), JHS 1 (t), JHS 2 (t) . . . , JHS n (t) were introduced and defined as: where ξ ∈ N was a natural number. The approximating candidate functions were then expressed as: where H m,i ∈ R were numeric coefficients related to the summation index m ∈ N. These coefficients, which were now not related to the amount of experimental information, were condensed into the H m vector. In this work, for the sake of simplicity, only the functions related to values of m ranging from 0 to 5 were explored (see Appendix A about how to build these functions).
More specifically, all the new functions in Equation (10) were used to numerically reproduce (nonlinear least squares method, Levemberg-Marquardt algorithm with all starting guess set to be zero) the evolution of the scar capsule thickness over time in MATLAB (MathWorks, Inc., Portola Valley, CA, USA), through the optimization of the numeric coefficients of the H m vector.
Furthermore, to quantify the intensity of the Runge phenomenon for each approximating function, the following dimensionless metric was used: where f was each of the previous approximating function (i.e., g 9 (t), g 10 (t), g 11 (t), g 12 (t), ) and max{s(t)} and min{s(t)} were the maximum and the minimum values of the experimental thickness. Finally, to figure out the overall suitability of each candidate function, the following dimensionless metric was introduced: where k = 10 4 [µm] was a dimensional constant, R 2 was the square of the multiple correlation coefficient, while SSE was the value of the sum of squared errors and d f = was the difference between the maximum and minimum values of the approximating function over time.

Results
The reconstruction of a continuous evolution of the scar tissue thickness from a small scattered set of experimental data through Equation (1) involved the use of the Vandermonde matrix in Equation (2). However, this matrix was badly conditioned (Figure 1a) and resulted in a polynomial function (Equation (3)) which was fully characterized by the elements of the a vectors (see Table 1). This polynomial was able to exactly reproduce experimental detections (i.e., R 2 = 1, SSE = 0), even if its evolution greatly differed from experimental points for inter-interval time values. In particular, the maximum difference between predicted values and expected ones was about 1.62 × 10 8 at t = 4839 h (Figure 1b).
Trying to improve this approach, the experimental times were normalized over the maximum time, then all the normalized experimental times ranged between 0 and 1. In this second case, the Vandermonde matrix was better conditioned (Figure 1c), but the function was still affected by the same numerical error in predictions (the maximum error was the same) even if scaled on abscissas (in this case the normalized time was 0.897), as shown in Figure 1d. Since the classic Vandermonde approach failed, a second approach was implemented to find a continuous function reproducing the scar tissue dynamics. First, Equation (4) was used to obtain a polynomial approximating function g 9 (t) (Figure 2a), which, by definition, was able to exactly reproduce all experimental points (R 2 = 1, SSE = 0). This function, which was expressed through Equation (5), was fully characterized by the x vector whose coefficients are equals to those of the a vector resulting from the previous approach (see Table 1). In particular, since the Runge phenomenon arose, at t = 4838 h the predicted thickness of the scar capsule was s(4838) = −1.62 × 10 8 µm.
To improve this standard approach, some further time points were added through the use of Equations (6) and (7). In particular, as a result of the first iteration, the time point t = 4032 h was added and the predicted thickness was s(4032) = 47.039 µm. By adding this time point to the experimental data set (i.e., new total number of detections n = 9 + 1 = 10), the resulting approximating function g 10 (t) (Figure 2b) was characterized by the x vector (see Table 2, column 1). In this case, at t = 5016 h, the predicted thickness was s(5016) = 3.656 × 10 7 µm because of the arising numerical oscillations.
Again, by iterating this process, an additional time point was created at t = 4704, while the predicted thickness was s(4704) = 47.025 µm. The coefficients of the approximating function g 11 (t) are shown in Table 2, column 2, while the maximum thickness value was s(5167) = −3.134 × 10 6 µm at t = 5167 due to the influence of the Runge phenomenon. Finally, by iterating the previous procedure once more, an additional time point was created at t = 5040 h, the predicted thickness was s(5040) = 47.018 µm and the numerical coefficients characterizing the candidate function g 12 (t) are shown in Table 2, column 3. Moreover, in this case, the Runge phenomenon arose leading to a prediction of s(3397) = −2.53 × 10 5 µm at t = 3397 h. It is worthy to note that for all the previous candidate functions g 10 (t), g 11 (t), g 12 (t), by definition, standard metrics values were R 2 = 1, SSE = 0.  In order to provide a novel way to reconstruct a continuous evolution from the experimental data set, Equation (10) was tested for m = 0, . . . , 5. In all cases, the elements of each H m vector are shown in Figure 3 (lines 0-5). In particular, for m = 0, the function was constant (H 0,0 = 22.51) and resulted in nonnegligible errors (i.e., SSE = 4181, see Figure 4a). For m = 1, instead, the function w 1 (t) was able to follow the experimental evolution and resulted in R 2 = 0.866, SSE = 347.3 (Figure 4b). Similarly, for m = 2, the function w 2 (t) led to R 2 = 0.9195, SSE = 208.9 (Figure 4c), while for m = 3 and then for the candidate function w 3 (t), the values of standard metrics were R 2 = 0.9197, SSE = 208.4 ( Figure 4d). Finally, for m = 4 and m = 5, the functions w 4 (t) and w 5 (t) resulted in  To assess the ability of each candidate function to reproduce the experimental set of data without large numerical oscillations, the intensity of the Runge phenomenon was quantified. More specifically, the performance of q 9 (t), achieved through the Vandermonde approach, was measured through the quantity ln(Λ q 9 (t) ) which resulted in 2.7188. Similarly, for g 9 (t), g 10 (t), g 11 (t), g 12 (t), obtained through the Lagrange approach, the values of the natural logarithm of the metric in Equation (11)  On the contrary, for the novel function achieved through Equation (10), the values of ln(Λ f ) were −1.3990, −3.2687, −3.6764, −3.8140 and −2.2977 for w 1 (t), w 2 (t), w 3 (t), w 4 (t) and w 5 (t), respectively, while it was, by definition, lim d w 0 (t) →0 ln(Λ w 0 (t) ) → −∞ for w 0 (t) (Figure 5a).
On the contrary, for the novel function achieved through Equation (10), the I metric resulted in 1.6722, 23.5523, 35.9416, 42.7525 and 8.8138 for w 1 (t), w 2 (t), w 3 (t), w 4 (t), w 5 (t), respectively, while it was lim d w 0 (t) →0 I w 0 (t) → 0 − for w 0 (t) (Figure 5b).  Figure 5. (a) The natural logarithm of the Λ f index was used to account for the intensity of the Runge phenomenon since it is able to graphically show whether approximating functions were affected by the Runge phenomenon. The first approximation obtained through the Vandermonde approach resulted in a positive value of ln(Λ g 9 (t) ). Similarly, for g 9 (t), · · · , g 12 (t), this logarithm was positive since the Runge phenomenon affected predictions. On the contrary, for w 1 (t), · · · , w 5 (t), this logarithm resulted in negative values and the Runge phenomenon did not affect numerical predictions of inter-interval values. For the function w 0 (t), the index ln(Λ w 0 (t) ) increased without bound towards −∞ (not shown within the plot). (b) Overall suitability: Considering both the ability of each candidate function to reproduce experimental values and the intensity of the Runge phenomenon affecting numerical predictions, the overall suitability shows that only the w n (t) family was able to reproduce experiments in a suitable way. Among the w n (t) functions, the w 4 (t) was also the most suitable to reproduce experimental data for inter-interval times.

Discussion
In computational medical physics, the evolution of a complex phenomenon should be often reproduced without clear knowledge of all physical details. The use of the polynomial interpolation based on the Vandermonde matrix is one of the most straightforward methods to obtain an interpolating polynomial from a set of experimental data points. Nevertheless, this standard approach failed to reproduce the evolution of the encapsulating tissue thickness around polymeric implants within the peripheral nerves. Indeed, the Vandermonde matrix in Equation (2) was badly conditioned since the time interval between the first and the last experimental observation was quite large. In this case, several elements were computed as zero when Equation (1) was solved because of the huge difference among matrix elements in Equation (2).
A first relevant question was related to the nature of the incorrect representation of the encapsulating tissue thickness evolution. Indeed, it was likely affected by the previous numerical effect (i.e., the Vandermonde matrix badly conditioned). To explore this point, all times were normalized with respect to the maximum one. In this way, the numeric interval between the first and the last experimental detection was scaled and the Vandermonde matrix in Equation (2) was better conditioned. However, the interpolant polynomials were still affected by the Runge phenomenon, as shown in Figure 1c,d. A first consequence was then that the lack of precision in inter-interval predictions was intrinsically due to the small size of the data set, together with the scattered nature of the experimental detections.
Moreover, a second standard approach based on the Lagrange approximating polynomial was used. Indeed, this method is, in general, more suitable in the case of a small set of experimental data. Unfortunately, Equation (4) still led to the same approximating polynomials as the previous approach. Moreover, in this case, the approximating function not only oscillated for inter-interval times, but also provided nonphysical results as negative values for the thickness.
Therefore, the initial set of experimental data was increased through supplementary data points. They were calculated through the iterative applications of Equations (6) and (7) to the last points of the time series and added to the initial data set since a final steady state was guessed. This procedure was able to decrease the amplitude of numerical oscillations from about 1.5 × 10 8 to about 2 × 10 5 (see Figure 2b-d). However, all interpolant polynomials, which were connected through the Neville's formula in Equation (8), were still affected by the Runge phenomenon.
However, the proposed strategy involving the increase of the initial set of points was not exactly designed to solve the initial problem. Indeed, just the whole set of experimental points was the only source of information. More specifically, extra points added spurious information about the dynamics of the scar tissue evolution. Nevertheless, this strategy was used to show that the increase of the cardinality of the set of points was also not able to eliminate numerical oscillations. In other words, this strategy was only proof of the concept to show the small influence of these extra points. In addition, since the oscillations of the interpolant polynomials were huge, the spurious information about the thickness of the extra points (see Equation (7)) also became negligible with respect to the oscillation amplitude. Moreover, other more effective procedures involving the choice and/or the addition of extra points to lower the Runge phenomenon (e.g., Chebyshev nodes [61,62] on the time axis) seemed unusable since no explicit formulas for the scar tissue dynamics were available in advance. Similarly, the use of splines (e.g., cubic splines [63,64]) seemed unable to provide an explicit expression to link all the points of the experimental set but only a piecewise function linking subsets of points. After the failure of the previously discussed standard approaches, a novel framework involving the use of Equation (10) was proposed.
More specifically, this novel set of approximating functions was able to disentangle the accuracy of the numerical predictions from the amount of experimental information and improve the approximation of the capsule thickness evolution over time. In particular, all the functions with m ranging from 1 to 5 closely approximated the experimental thickness evolution since the standard statistics R 2 ranged from 0.8 to about 0.9 (a particular unsuitable case was found for m = 0 since w 0 (t) = 22.51 identically).
Nevertheless, the intensity of the Runge phenomenon was not described by standard metrics, then it was accounted through Equation (11). Indeed, this numerical phenomenon resulted in large oscillations between experimental data points, thus, even if standard metrics were maximized (i.e., R 2 = 1, SSE = 0), inter-interval predictions were not physically meaningful. As shown in Figure 5a, the natural logarithm of this metric was also able to visualize the onset of the Runge phenomenon. More specifically, a positive finite value of ln(Λ f ) was related to the onset of the Runge phenomenon while a finite negative value was related to the absence of this numerical effect. The function w 0 (t), for which ln(Λ w 0 (t) ) → −∞, was obviously unsuitable to reproduce experiments since it was constant.
Finally, the overall performances of each candidate function were compared through the use of Equation (12). In particular, this metric was directly proportional to the ability of each candidate function to explain the variance of experimental data while it was inversely proportional to the difference between the maximum and the minimum values of the approximating function. Similarly, it was inversely proportional to the sum of the square errors increased of 1 (to avoid inconsistent values when SSE = 0, e.g., for Lagrange polynomials), as well as to the Λ f metric, accounting for the onset of the Runge phenomenon. Through the use of Equation (12), the best candidate function was found for the family described in Equation (10) with m = 4.

Conclusions
This work provides an original procedure to reproduce the evolution of the mean thickness of the scar tissue encapsulating an implanted polymeric intraneural dummy electrode over time through a linear combination of novel functions. This approach allowed the experimental data to be reproduced in an effective way (R 2 = 0.92), avoiding the onset of the Runge phenomenon, keeping the residuals small (SSE = 202) and providing a good stability for inter-interval predictions. Furthermore, all the provided results were physically consistent (i.e., positive thickness) and comparable to experimental detections. In addition, the approximating function was smooth. Even if other approaches have been proposed to find approximating functions able to reconstruct the evolution of the scar capsule thickness over time [65], their performance is clearly overcome by our novel framework (i.e., R 2 = 0.86 vs. R 2 = 0.92 here). Moreover, since this novel approach is general, it could be used to numerically reproduce the scar tissue formation in different districts of the body where biomaterials are in contraction with biological tissue, even if, in this work, it was tested within in vivo peripheral nerves. Funding: This work was partially supported by the NeuHeart project (a neuroprosthesis to restore the vagal-cardiac closed-loop connection after heart transplantation-GA no. 824071). Publication funds were provided by Pier Nicola Sergi.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Conflicts of Interest:
The auhtors declare no conflict of interest.