## 3 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\times 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\times 10^{8}\) um. Figure 1: Vandermonde-based interpolation (black squares = mean experimental values, solid lines = theoretical predictions): (**a**) The huge difference among the element 9,9 and the other ones resulted in an ill-conditioned matrix. (**b**) The interpolating polynomial function was affected by the Runge phenomenon, resulting in nonphysical predictions for inter-interval values (negative values). (**c**) The normalization of times over the maximum experimental time lead to a better conditioned matrix. (**d**) Moreover, in the case of time normalization, the performances of the interpolating polynomials were unchanged. 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 \(\upmu\)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 \(\mathbf{x}\) vector (see Table 2, column 1). In this case, at \(t=5016\) h, the predicted thickness was \(s(5016)=3.656\times 10^{7}\)\(\upmu\)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\)\(\upmu\)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\times 10^{6}\)\(\upmu\)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\)\(\upmu\)m and the numerical coefficients characterizing the candidate function \(g_{12}(t)\) are shown in Table 2, column 3.