Novel In Silico Strategies to Model the In Vivo Nerve Scarring Around Implanted Parylene C Devices

Abstract

The implantation of materials into in vivo peripheral nerves triggers the production of scar tissue. A scar capsule progressively incorporates foreign bodies, which become insulated from the surrounding environment. This phenomenon is particularly detrimental in the case of electrical active sites enveloped within scar sheets, since the loss of contact with axons highly decreases the effectiveness of neural interfaces. As a consequence, the in silico modelling of scar capsule evolution may lead to improvements in the design of intraneural structures and enhancing their reliability over time. In this work, a novel theoretical framework is proposed to model the evolution of capsule thickness over time together with an improved optimisation procedure able to avoid apparently suitable choices resulting from standard procedures. This framework provides a fast, simple, and accurate modelling of experimental data ($R^2=0.97$), definitely improving on previous approaches.

1. Introduction

Peripheral nerves provide the natural connection between the periphery of the body and the brain [1]. They are extremely complex structures, which, nonetheless, can be connected to external high-tech devices to record neural signals and evoke neural activity. These devices, which are commonly referred as neural interfaces, have been designed to be implanted around or directly inserted within the peripheral nerves [2,3,4,5], optionally through microneedles [6,7,8]. Similarly, other kinds of electrodes (i.e., regenerative neural interfaces) were investigated [2,9,10] and found to benefit from the ability of axons to regenerate [11,12,13,14]. In particular, intraneural interfaces, despite their relative invasiveness, were able to provide selective recording and stimulation.

However, the foreign body reaction (FBR) normally affects the interaction between living tissues and synthetic materials [15,16]. In this case, the FBR influences the interaction between peripheral nerve tissue [1,17,18,19,20,21,22,23] and the intraneural interface [24,25,26], leading to the formation of a tissue capsule [27,28] around the implanted device that hampers its functionality. Several works investigated the FBR and the main phases were described: first, blood and plasma proteins (e.g., fibrinogen, fibronectin, albumin, and antibodies) are absorbed to the implant surface [29,30]; then, the adsorbed protein layer promotes monocyte and leukocyte extravasation, attraction, and adhesion to the surface, together with the coagulation cascade [31,32,33]. Furthermore, a matrix of fibrin was found to attract leukocytes and macrophages around the surface of the implanted material through the chemical attraction of chemokines [32,34]. Again, the extravasated monocytes differentiate to macrophages and fuse together to form multinucleated foreign body giant cells (FBGCs), which start releasing further inflammatory cytokines, boosting the inflammatory response through a positive feedback mechanism [35,36]. Thus, the myofibroblasts trigger the secretion of extracellular matrix (ECM) components [16,27,37]. Finally, the scar tissue capsule was found to be impermeable to the nonspecific immune system and to many chemicals [15,16,27].

Currently, a reliable prediction of the scar capsule thickness around an inserted material is difficult because of the complexity of the FBR [15,16,27,28,29,30,31,32,33,34,35,36,37,38] and the unclear physical and mechanical interactions with peripheral nerves [8,39,40,41,42,43,44,45,46,47,48,49]. Previous literature investigated the FBR evolution triggered by intraneural electrodes longitudinally implanted (e.g., LIFE electrodes [50,51]) and studied the cellular composition and the thickness of the fibrotic capsule up to 8 months [52,53]. These studies found that the increase in the thickness of the scar capsule was directly related to the increase in the electrical impedance, which resulted in the decrease in the signal-to-noise ratio [54], together with a drop in recording and stimulation capabilities [36,55].

Again, the quantitative prediction of the scar tissue thickness emerged as crucial knowledge to improve both the design and the implantation techniques of intraneural interfaces, even if several interrelated FBR processes are partially hidden and their biophysical rules are not known in detail [15,16,27,28,29,30,31,32,33,34,35,36,37,38]. Therefore, there is a pressing issue to investigate what kind of computational strategies could be able to model scar tissue evolution, also for small experimental data sets with close experimental time points, which are prone to be ill conditioned. In addition, the high complexity and the huge number of biophysical and biochemical processes may result in extremely complex and time-consuming computational models.

As a consequence, in this work a novel phenomenological model is proposed to achieve a physically compatible evolution of the mean capsule thickness together with time-saving simulations able to start from a set of scattered experimental data with low cardinality [56]. Although the proposed phenomenological approach was previously successfully used [57,58], here the previous numerical performances are clearly improved upon, while the computational complexity is kept low. More specifically, this work is organised as follows: first, the novel model is presented, and then an optimisation procedure, involving the use of standard statistics, is performed as a benchmark. The suitable candidate functions, selected in the previous analysis, are, then, validated with respect to their ability to model the experimental data set with standard levels of statistical significance. Again, they are further tested with respect to their physical consistency (in order to avoid nonphysical behaviours). Thus, a novel metric is provided to identify the best form of the candidate function. A sensitivity analysis of this function is performed to investigate how the global ability to model experimental data is dependent on selected parameters, while a final optimisation of all parameters of the best functional form is performed.

2. Methods

Parylene C double devices (20 mm long, 200 $\mathsf{\mu }$m wide, and 10 $\mathsf{\mu }$m thick) were longitudinally implanted within the tibial branch of the sciatic nerve of female Sprague Dawley rats (250–300 g, n = 6–8 for each group) [56] and resulted in the formation of a scar tissue capsule. The thickness of this capsule formed around these implanted Parylene C devices was examined at 0, 24, 48, 96, 336, 672, 1344, 2688, and 5376 h post-implant by light microscopy in $4x$ images and measured as the distance between each side of the device and the closest myelinated axon though ImageJ software (version v1.51, National Institute of Health, USA).

The evolution of the scar tissue thickness over time was assumed to be modelled through an unknown continuum function $s\left(t\right):A\to B$, where $A,B\subseteq R^{+}$. However, here, the subset $S=\{s\left(t_1\right),\dots ,s\left(t_9\right)\}$ of mean experimental values (with standard deviation) coming from the whole experimental set [56] was used. Let $P_i(\overline{\alpha },t)$ be a polynomial function defined as follows:

$$\begin{array}{c}\hfill P_i(\overline{\alpha },t)=\{\sum _{i=0}^n{\alpha }_i{t}^i|{\alpha }_0,\dots ,{\alpha }_n\in R\}\end{array}$$

(1)

while two supplementary functions are defined as follows:

$$\begin{array}{c}\hfill g\left(t\right)=abs\left(t\right)\end{array}$$

(2)

$$\begin{array}{c}\hfill h\left(t\right)=ln\left(t\right)\end{array}$$

(3)

where t is the time, while $abs\left(t\right)$ and $ln\left(t\right)$ are the absolute value and the natural logarithm of t, respectively. Therefore, the operator M is defined as follows:

$$M[P_i(\overline{\alpha },t),F_j(\overline{\beta },t),\gamma ]=P_i(\overline{\alpha },t)+\gamma F_j(\overline{\beta },t)$$

(4)

where $\gamma \in \mathbb{R}$ and i and j are the maximum degree of the polynomial functions. The vectors $\overline{\alpha }={[{\alpha }_0,\dots ,{\alpha }_i]}^T$ and $\overline{\beta }={[{\beta }_0,\dots ,{\beta }_j]}^T$ include the coefficients of the polynomial functions. In addition,

$$\begin{array}{c}\hfill F_j(\overline{\beta },t)=(h\circ g\circ P_j)(\overline{\beta },t)\end{array}$$

(5)

where “∘” is the composition operator among functions. As a consequence, the evolution over time of the scar tissue thickness was modelled through the approximated function:

$$\begin{array}{c}\hfill q(t,i,j,\overline{\alpha },\overline{\beta },\gamma )=M[P_i(\overline{\alpha },t),F_j(\overline{\beta },t),\gamma ]\end{array}$$

(6)

Because of the polynomial nature of $P_i(\overline{\alpha },t)$ and the cardinality of the experimental data set ($n=9$), the minimum degree of polynomials resulting in Runge instability was $i=5$ [57]. Thus, in this work, the attention was narrowed to candidate functions with $i=1,\dots ,5$ and $j=1,\dots ,5$, while, in general, $i\ne j$ and $\overline{\alpha }\ne \overline{\beta }$. However, to reduce the number of elements of vectors $\overline{\alpha },\overline{\beta }$ the candidate functions were chosen in order to satisfy the further constraints ${\alpha }_i={\beta }_i$ for $i\le j$. To test the ability of Equation (6) to model the experimental set of data, the standard statistics $RMSE(i,j),R^2(i,j)$ were computed for each set of best $\overline{\alpha },\overline{\beta }$, resulting from a standard optimisation procedure (nonlinear least-squares method, least absolute residuals, Levemberg–Marquardt algorithm) in Matlab (Mathworks Inc. Academic licence). In particular, each element of a 5$\times$5 matrix W was defined as follows:

$$\begin{array}{c}\hfill w(i,j)=lo{g}_{10}\left[RMSE(i,j)\right]\end{array}$$

(7)

where each $RMSE(i,j)$ was the value of the $RMSE$ (Round Mean Square Error) related to the fitting performances of the candidate function $q(t,i,j,\overline{\alpha },\overline{\beta })$, defined for the best set of parameters $\overline{\alpha }$,$\overline{\beta }$. Therefore, each element of a twin matrix $LRR$ was defined as follows:

$$\begin{array}{c}\hfill lrr(i,j)=lo{g}_{10}\left[R^2(i,j)\right]\end{array}$$

(8)

while the primary matrix $RR$ was defined through

$$\begin{array}{c}\hfill rr(i,j)=R^2(i,j)\end{array}$$

(9)

where each $R^2(i,j)$ was the value of the standard $R^2$ (coefficient of determination) related to the fitting performances of the candidate function $q(t,i,j,\overline{\alpha },\overline{\beta })$, defined for the best set of parameters $\overline{\alpha }$,$\overline{\beta }$.

All combinations of i and j resulting in small values of $w(i,j)$ were related to potentially suitable candidate functions, while those resulting in high values of $w(i,j)$ were related to unsuitable candidate functions. The W matrix was written as follows:

$$\begin{gathered} W=\left(\begin{array}{ccccc}{W}_{A1} \\ {W}_{B1}& W_{A2} \\ {W}_{B2}& W_{A3} \\ {W}_{B3}& W_{A4} \\ {W}_{B4}& W_{A5} \\ {W}_{B5}\end{array}\right) \end{gathered}$$

(10)

where $W_{A\lambda }$, $\lambda =1,\dots ,5$, were block matrices of potentially suitable combinations, while $W_{B_{\lambda }}$ were block matrices of unsuitable combinations of indexes $i,j$. Similarly, for the twin matrix, defined through Equation (8),

$$LRR=\left(\begin{array}{ccccc}LR{R}_{A1}& LR{R}_{A2}& LR{R}_{A3}& LR{R}_{A4}& LR{R}_{A5}\\ LR{R}_{B1}& LR{R}_{B2}& LR{R}_{B3}& LR{R}_{B4}& LR{R}_{B5}\end{array}\right)$$

(11)

The same structure of the twin matrix was finally applied to the $RR$ matrix, defined through Equation (9), which was written as follows:

$$RR=\left(\begin{array}{ccccc}R{R}_{A1}& R{R}_{A2}& R{R}_{A3}& R{R}_{A4}& R{R}_{A5}\\ R{R}_{B1}& R{R}_{B2}& R{R}_{B3}& R{R}_{B4}& R{R}_{B5}\end{array}\right)$$

(12)

where $R{R}_{A\lambda }$, $\lambda =1,\dots ,5$, were block matrices related to potentially suitable combinations, while $R{R}_{B\lambda }$ were block matrices of unsuitable combinations of indexes $i,j$. A comparison between experimental data and in silico predictions coming from each candidate function was used to validate each potentially suitable candidate function. The distribution of in silico values coming from each candidate function was tested through the Shapiro–Francia normality test with a confidence level of $0.05$, while the statistical significance of the difference between experimental data and in silico predictions was tested through an unpaired Student t-test with both $0.05$ and $0.01$ confidence levels.

Nevertheless, since the use of Equation (6) to model the scar tissue evolution was novel, the global effectiveness of the previous standard procedure was not totally expected. Therefore, further tests were performed to investigate whether nonphysical effects were hidden in potentially suitable candidate functions. With this aim, each single block was analysed starting from $R{R}_{A5}$ to end with the block $R{R}_{A1}$. More specifically, $n_d(i,j)=card\left\{N_d(i,j)\right\}$ was the cardinality of the set $N_d(i,j)$, which for each candidate function accounted for the discontinuity points

$$N_d(i,j)=\left\{t_{\eta }\right|li{m}_{t\to t_{\eta }^{-}}q(t,i,j,\overline{\alpha },\overline{\beta },\gamma )\ne li{m}_{t\to t_{\eta }^{+}}q(t,i,j,\overline{\alpha },\overline{\beta },\gamma )\}$$

(13)

or accounted for the vertical asymptotes

$$N_d(i,j)=\left\{t_{\eta }\right|li{m}_{t\to t_i^{\pm ,\mp }}q(t,i,j,\overline{\alpha },\overline{\beta },\gamma )=\pm \infty \}$$

(14)

with $\eta \in N$. Furthermore, $n_z(i,j)=card\left\{N_z(i,j)\right\}$ was the cardinality of the set $N_z(i,j)$, which accounted for the stationary points and it was defined as:

$$N_z(i,j)=\left\{t_{\eta }\right|\dot{q}(t,i,j,\overline{\alpha },\overline{\beta },\gamma )=0\}$$

(15)

Starting from the previous analysis, each element of a novel matrix $\mathsf{\Xi }(RR,n_d,n_z)$ was defined for each combination of i and j as follows:

$$\xi [i,j,rr(i,j),n_d(i,j),n_z(i,j)]=rr(i,j)[1-n_d(i,j)][1-n_z(i,j)/10]$$

(16)

Each element of the matrix $\xi [i,j,rr(i,j),n_d(i,j),n_z(i,j)]$ was a novel metric, able to evaluate at the same time both the standard performance and the presence of numerical singularities for each trial function. This matrix was used to identify the best functional form of the candidate function, resulting from the best combination of $i,j$ indexes $i^{*}$ and $j^{*}$. In addition, to further explore this functional form, the changes in shape were analysed for numerical changes in parameters $\overline{\alpha },\overline{\beta },\gamma$, and a sensitivity analysis was performed. More specifically, the increments $\Delta \overline{\alpha },\Delta \overline{\beta },\Delta \gamma$ were used, where, in general, $\Delta \overline{\alpha }={[\Delta {\alpha }_0,\dots ,\Delta {\alpha }_i]}^T$, $\Delta \overline{\beta }={[\Delta {\beta }_0,\dots ,\Delta {\beta }_j]}^T$, $\Delta {\alpha }_i=\Delta c$, $\Delta {\beta }_j=\Delta c$, $\Delta \gamma =\Delta c$, for each relevant value of the $i,j$ indexes at a time, and $\Delta c=0,1,5,10,50,100,500,1000$, were the different increments. The sensitivity of the functional form with respect to each parameter was, then, defined as follows [59,60]:

$$S\left[q(t,i^{*},j^{*},\overline{\alpha }+\Delta \overline{\alpha },\overline{\beta },\gamma )\right]={\displaystyle \frac{q(t,i^{*},j^{*},\overline{\alpha }+\Delta \overline{\alpha },\overline{\beta },\gamma )-q(t,i^{*},j^{*},\overline{1},\overline{1},1)}{q(t,i^{*},j^{*},\overline{\alpha }+\Delta \overline{\alpha },\overline{\beta },\gamma )}}$$

(17)

$$S\left[q(t,i^{*},j^{*},\overline{\alpha },\overline{\beta }+\Delta \overline{\beta },\gamma )\right]={\displaystyle \frac{q(t,i^{*},j^{*},\overline{\alpha },\overline{\beta }+\Delta \overline{\beta },\gamma )-q(t,i^{*},j^{*},\overline{1},\overline{1},1)}{q(t,i^{*},j^{*},\overline{\alpha },\overline{\beta }+\Delta \overline{\beta },\gamma )}}$$

(18)

$$S[q(t,i^{*},j^{*},\overline{\alpha },\overline{\beta },\gamma +\Delta \gamma ]={\displaystyle \frac{q(t,i^{*},j^{*},\overline{\alpha },\overline{\beta },\gamma +\Delta \gamma )-q(t,i^{*},j^{*},\overline{1},\overline{1},1)}{q(t,i^{*},j^{*},\overline{\alpha },\overline{\beta },\gamma +\Delta \gamma )}}$$

(19)

Finally, numeric instability effects due to the change in the $\gamma$ parameters were studied with reference to the variation in the modelling ability of the best candidate function. The best value of the $\gamma$ parameter was found through the optimisation of the standard $R^2$ statistic.

3. Results

3.1. A Standard Approach to Discriminate Suitable Candidate Functions

A first discrimination between suitable and not suitable trial functions was performed through the W matrix. More specifically, for the family of functions with $i=1$ the values of $w(1,j)$ (with $1\le j\le 5$) ranged from $0.836$ to $1.475$, while for functions with $i=2$ the values $w(2,j)$ varied from $0.805$ to $3.029$. Similarly, for $i=3$ the $w(3,j)$ values varied from $0.837$ to $3.059$, while for functions with $i=4$ the range of variation in $w(4,j)$ was between $0.686$ and $3.068$. Finally, for candidate functions with $i=5$ the values $w(5,j)$ varied between $0.943$ and $3.205$. Moreover, the main 5$\times$5 W matrix was divided into two different blocks, as shown in Figure 1a. The upper block matrix was formed by five submatrices $W_{A1}$ (4$\times$1), $W_{A2}$ (3$\times$1), $W_{A3}$ (3$\times$1), $W_{A4}$ (3$\times$1), and $W_{A5}$ (3$\times$1), and it was related to combinations of i and j indexes resulting in candidate functions potentially suitable to model the experimental set of data (Figure 1b). On the contrary, the lower block matrix was formed by five submatrices $W_{B1}$(1$\times$1), $W_{B2}$ (2$\times$1), $W_{B3}$ (2$\times$1), $W_{B4}$ (2$\times$1), and $W_{B5}$ (2$\times$1), and it was related to the combinations of i and j indexes, which resulted in candidate functions unsuitable to model the experimental set of data (Figure 1c).

Equivalently, within the LRR matrix (decimal logarithm of the coefficient of determination) for the family of functions with $i=1$ the values of $lrr(1,j)$ (with $1\le j\le 5$) ranged from $-0.029$ to $-0.056$, while for functions with $i=2$ the values $lrr(2,j)$ varied from $-0.026$ to $-0.041$. Furthermore, for $i=3$ the $lrr(3,j)$ values varied from $-0.021$ to $-0.046$, while for functions with $i=4$ the range of variation in $lrr(4,j)$ was between $-0.027$ and $3.297$. Finally, for candidate functions with $i=5$ the values $lrr(5,j)$ varied between $-0.503$ and $0.417$.

Again, the main 5$\times$5 $LRR$ matrix was divided into two different blocks, as shown in Figure 1d. The upper block matrix was formed by five submatrices $LR{R}_{A1}$ (4$\times$1), $LR{R}_{A2}$ (3$\times$1), $LR{R}_{A3}$ (3$\times$1), $LR{R}_{A4}$ (3$\times$1), and $LR{R}_{A5}$ (3$\times$1), and it was related to combinations of i and j indexes resulting in candidate functions apparently suitable to model the experimental set of data (Figure 1e). The lower block matrix was formed by the submatrices $LR{R}_{B1}$(1$\times$1), $LR{R}_{B2}$ (2$\times$1), $LR{R}_{B3}$ (2$\times$1), $LR{R}_{B4}$ (2$\times$1), and $LR{R}_{B5}$ (2$\times$1), which were related to the combinations of i and j indexes resulting in functions not suitable to model the experimental set of data (Figure 1f).

The same division in submatrices was applied to the matrix $RR$, resulting in the following blocks, $R{R}_{A1}$(4$\times$1), $R{R}_{A2}$ (3$\times$1), $R{R}_{A3}$ (3$\times$1), $R{R}_{A4}$ (3$\times$1), and $R{R}_{A5}$ (3$\times$1), which were related to suitable combinations of i and j indexes. Correspondingly, the lower block matrix was formed by $R{R}_{B1}$ (1$\times$1), $R{R}_{B2}$ (2$\times$1), $R{R}_{B3}$ (2$\times$1), $R{R}_{B4}$ (2$\times$1), and $R{R}_{B5}$ (2$\times$1). The latter blocks (B type) were excluded from further investigations, while the former ones were explored in more depth, starting from the $R{R}_{Ai}$ blocks with higher values.

3.2. Standard and Not Standard Analysis of Different Families of Candidate Functions

The in silico predictions resulting from all suitable combinations were compared to experimental data to validate them. In addition, the form of their distribution was tested for normality and compared to experimental data to quantify eventually significant differences. Then, the continuum evolution of each candidate function was explored to reveal eventual nonphysical behaviours.

In particular, the $RR$ value referring to the candidate function with $i=1$ and $j=5$ was quite elevated ($R{R}_{A5}\left(1\right)=0.9343$), while the residuals between experimental and predicted values oscillated between $-5.385$ and $10.641$ $\mathsf{\mu }$m (Table S1 and Figure S1). This function was able to show a good positive correlation between predictions and experimental data ($R^2=0.94$, Figure 2a). The Shapiro–Francia normality test resulted in $p=0.94$ for in silico values, while in $p=0.48$ for experimental data (Figure 2b). Furthermore, an unpaired Student t-test between experimental data and in silico predictions resulted in $p=0.86$ (Figure 2c). Nevertheless, a complete plot of the function within the experimental time range showed two vertical asymptotes at $t=24.000$ and at $t=47.338$ h, together with three stationary points at $t=0$, $t=17.106$, and $t=40.171$ h (Figure 2d and Figure 2d inset). Then, recalling Equations (13)–(15), for this function, $n_d=2$ and $n_z=3$.

The $RR$ value associated with the candidate function with $i=2$ and $j=5$ was high ($R{R}_{A5}\left(2\right)=0.9406$), while the residuals between experimental and predicted values oscillated between $-5.963$ and $8.414$ $\mathsf{\mu }$m (Table S2 and Figure S2). This function provided a similar correlation ($R^2=0.94$) between experimental and predicted data (Figure 2e). The Shapiro–Francia normality test resulted in $p=0.57$ for in silico values (Figure 2f), while an unpaired Student t-test between experimental data and in silico predictions resulted in $p=0.99$ (Figure 2g). However, the complete plot showed two vertical asymptotes at $t=24.000$ and $t=47.338$ h, together with three stationary points at $t=0$, $t=17.106$, and $t=40.171$ h (Figure 2h and Figure 2h inset). Therefore, for this function $n_d=1$ and $n_z=1$.

Similarly, the $RR$ value pertaining to the candidate function deriving from the combination of $i=3$ and $j=5$ was elevated ($R{R}_{A5}\left(3\right)=0.952$), and the residuals between experimental and predicted values oscillated between $-4.171$ and $4.955$ $\mathsf{\mu }$m (Table S3 and Figure S3). This function showed a quite high correlation ($R^2=0.98$) between experimental and predicted data (Figure 2i). The Shapiro–Francia normality test resulted in $p=0.48$ for in silico values (Figure 2j), while an unpaired Student t-test between experimental data and in silico predictions resulted in $p=0.98$ (Figure 2k). A graphical representation of the function showed a vertical asymptote at $t=19.539$ h, together with a stationary point at $t=4218.948$ h (Figure 2l and Figure 2l inset). As a consequence, for this function $n_d=1$ and $n_z=1$. The numeric values of the $\overline{\alpha }$ and $\overline{\beta }$ elements are detailed in

Table 1. The numeric values of each element of $\overline{\alpha }$ and $\overline{\beta }$ vectors for the functions $q(t,1,5,\overline{\alpha },\overline{\beta },\gamma )$, $q(t,2,5,\overline{\alpha },\overline{\beta },\gamma ),$ and $q(t,3,5,\overline{\alpha },\overline{\beta },\gamma )$ are listed in columns 1, 2, and 3, respectively.

$\mathit{q}(\mathit{t},1,5,\overline{\mathit{\alpha }}$,$\overline{\mathit{\beta }}$,$\mathit{\gamma }$)	$\mathit{q}(\mathit{t},2,5,\overline{\mathit{\alpha }}$,$\overline{\mathit{\beta }}$,$\mathit{\gamma }$)	$\mathit{q}(\mathit{t},3,5,\overline{\mathit{\alpha }}$,$\overline{\mathit{\beta }}$,$\mathit{\gamma }$)
${\alpha }_0={\beta }_0=-1.433·{10}^{-03}$	${\alpha }_0={\beta }_0=-7.820·{10}^{-00}$	${\alpha }_0={\beta }_0=-3.611·{10}^{-07}$
${\alpha }_1={\beta }_1=1.021·{10}^{-01}$	${\alpha }_1={\beta }_1=4.718·{10}^{-02}$	${\alpha }_1={\beta }_1=-1.314·{10}^{-05}$
${\beta }_2=-1.619·{10}^{-00}$	${\alpha }_2={\beta }_2=4.964·{10}^{-02}$	${\alpha }_2={\beta }_2=1.825·{10}^{-09}$
${\beta }_3=-1.433·{10}^{-01}$	${\beta }_3=-1.765·{10}^{-06}$	${\alpha }_3={\beta }_3=-1.439·{10}^{-05}$
${\beta }_4=2.768·{10}^{-03}$	${\beta }_4=1.487·{10}^{-02}$	${\beta }_4=2.647·{10}^{-02}$
${\beta }_5=-7.636·{10}^{-01}$	${\beta }_5=-6.130·{10}^{-01}$	${\beta }_5=2.432·{10}^{-00}$

for the functions $q(t,1,5,\overline{\alpha },\overline{\beta },\gamma ),$ $q(t,2,5,\overline{\alpha },\overline{\beta },\gamma ),$ and $q(t,3,5,\overline{\alpha },\overline{\beta },\gamma )$.

The family of candidate functions with $j=4$ was then tested. The $RR$ value related to the candidate function with $i=1$ and $j=4$ was high ($R{R}_{A4}\left(1\right)=0.9272$), while the residuals between experimental and predicted values oscillated between $-5.519$ and $5.457$ $\mathsf{\mu }$m at the experimental times (Table S4 and Figure S4). This function was able to show a good positive correlation between predictions and experimental data ($R^2=0.97$, Figure 3a), and the Shapiro–Francia normality test resulted in $p=0.69$ for in silico values (Figure 3b). Furthermore, an unpaired Student t-test between experimental data and in silico predictions resulted in $p=0.96$ (Figure 3c). Nevertheless, the complete plot of the function in the experimental range showed three vertical asymptotes at $t=10.872$, $t=30.972$, and $t=659.550$ h, together with two stationary points at $t=22.453$ and $t=491.255$ h (Figure 3d and Figure 3d inset). Then, for this function $n_d=3$ and $n_z=2$.

The $RR$ value referring to the candidate function with $i=2$ and $j=4$ was also elevated ($R{R}_{A4}\left(2\right)=0.95$), while while the residuals between experimental and predicted values oscillated between $-4.564$ and $9.647$ $\mathsf{\mu }$m (Table S5 and Figure S5). This candidate showed a positive correlation ($R^2=0.95$) between experimental and predicted data (Figure 3e). The Shapiro–Francia normality test resulted in $p=0.36$ for in silico values (Figure 3f), while an unpaired Student t-test between experimental data and in silico predictions resulted in $p=0.94$ (Figure 3g). However, the complete plot showed a vertical asymptote at $t=23.995$ h, together with two stationary points at $t=18.867$ and $t=4023.990$ h (Figure 3h and Figure 3h inset). Therefore, for this function $n_d=1$ and $n_z=2$.

Equivalently, the $RR$ value pertaining to the candidate function with $i=3$ and $j=4$ was considerable ($R{R}_{A4}\left(3\right)=0.9538$), while the residuals between experimental and predicted values oscillated between $-4.679$ and $4.469$ $\mathsf{\mu }$m (Table S6 and Figure S6). This function resulted in a positive high correlation ($R^2=0.98$) between experimental and predicted data (Figure 3i), while the Shapiro–Francia normality test resulted in $p=0.50$ for in silico values (Figure 3j). Again, an unpaired Student t-test between experimental data and in silico predictions resulted in $p=0.99$ (Figure 3k). The complete plot of the function was able to underline a vertical asymptote at $t=19.627$ h, together with a stationary point at $t=4217.266$ h (Figure 3l and Figure 3l inset). As a consequence, for this function $n_d=1$ and $n_z=1$. The numeric values of the $\overline{\alpha }$ and $\overline{\beta }$ elements are listed in

Table 2. The numeric values of the elements of the $\overline{\alpha }$ and $\overline{\beta }$ vectors for the functions $q(t,1,4,\overline{\alpha },\overline{\beta },\gamma ),$ $q(t,2,4,\overline{\alpha },\overline{\beta },\gamma ),$ and $q(t,3,4,\overline{\alpha },\overline{\beta },\gamma )$ are shown in columns 1, 2, and 3, respectively.

$\mathit{q}(\mathit{t},1,4,\overline{\mathit{\alpha }},\overline{\mathit{\beta }},\mathit{\gamma }$)	$\mathit{q}(\mathit{t},2,4,\overline{\mathit{\alpha }},\overline{\mathit{\beta }},\mathit{\gamma }$)	$\mathit{q}(\mathit{t},3,4,\overline{\mathit{\alpha }},\overline{\mathit{\beta }},\mathit{\gamma }$)
${\alpha }_0={\beta }_0=-1.433·{10}^{-03}$	${\alpha }_0={\beta }_0=-7.820·{10}^{-00}$	${\alpha }_0={\beta }_0=-3.611·{10}^{-07}$
${\alpha }_1={\beta }_1=1.021·{10}^{-01}$	${\alpha }_1={\beta }_1=4.718·{10}^{-02}$	${\alpha }_1={\beta }_1=-1.314·{10}^{-05}$
${\beta }_2=-1.619·{10}^{-00}$	${\alpha }_2={\beta }_2=4.964·{10}^{-02}$	${\alpha }_2={\beta }_2=1.825·{10}^{-09}$
${\beta }_3=-1.433·{10}^{-01}$	${\beta }_3=-1.765·{10}^{-06}$	${\alpha }_3={\beta }_3=-1.439·{10}^{-05}$
${\beta }_4=2.768·{10}^{-03}$	${\beta }_4=1.487·{10}^{-02}$	${\beta }_4=2.647·{10}^{-02}$
${\beta }_5=-7.636·{10}^{-01}$	${\beta }_5=-6.130·{10}^{-01}$	${\beta }_5=2.432·{10}^{-00}$

for the functions $q(t,1,4,\overline{\alpha },\overline{\beta },\gamma ),$$q(t,2,4,\overline{\alpha },\overline{\beta },\gamma ),$ and $q(t,3,4,\overline{\alpha },\overline{\beta },\gamma )$.

Moreover, the family of candidate functions with $j=3$ was tested. The $RR$ value associated with the candidate function with $i=1$ and $j=3$ was elevated ($R{R}_{A3}\left(1\right)=0.9271$), while the residuals between experimental and predicted values oscillated between $-4.248$ and $6.877$ $\mathsf{\mu }$m at the experimental times (Table S7 and Figure S7). This function was able to show a good positive correlation between predictions and experimental data ($R^2=0.94$, Figure 4a), and the Shapiro–Francia normality test resulted in $p=0.88$ for in silico values (Figure 4b). Furthermore, an unpaired Student t-test between experimental data and in silico predictions resulted in $p=0.82$ (Figure 4c). Nevertheless, the function had a vertical asymptote at $t=21.176$ h, together with two stationary points at $t=0$ and $t=12.134$ h (Figure 4d and Figure 4d inset). Then, for this function $n_d=1$ and $n_z=2$.

The $RR$ value related to the trial function with $i=2$ and $j=3$ was also high ($R{R}_{A3}\left(2\right)=0.92$), while the residuals between experimental and predicted values oscillated between $-10.530$ and $3.007$ $\mathsf{\mu }$m (Table S8 and Figure S8). This function was able to show a good positive correlation between predictions and experimental data ($R^2=0.94$, Figure 4e), and the Shapiro–Francia normality test resulted in $p=0.46$ for in silico values (Figure 4f). Furthermore, an unpaired Student t-test between experimental data and in silico predictions resulted in $p=0.85$ (Figure 4g). However, the complete plot showed a vertical asymptote at $t=18.708$ h, together with two stationary points at $t=0.207$ and $t=4000.265$ h (Figure 4h and Figure 4h inset). Therefore, for this function $n_d=1$ and $n_z=2$.

Correspondingly, the $RR$ value pertaining to the trial function with $i=3$ and $j=3$ was considerably lower ($R{R}_{A3}\left(3\right)=0.4345$), while the residuals between experimental and predicted values oscillated between $-2.250$ and $11.407$ $\mathsf{\mu }$m (Table S9 and Figure S9). This function was able to show a highly positive correlation between predictions and experimental data ($R^2=0.98$, Figure 4i), and the Shapiro–Francia normality test resulted in $p=0.29$ for in silico values (Figure 4j). Furthermore an unpaired Student t-test between experimental data and in silico predictions resulted in $p=0.77$ (Figure 4k). The complete plot of the function revealed a stationary point at $t=4747.488$ h (Figure 4l). As a consequence, for this function $n_d=0$ and $n_z=1$. The numeric values of the $\overline{\alpha }$ and $\overline{\beta }$ elements are detailed in

Table 3. The numeric values of the elements of $\overline{\alpha }$ and $\overline{\beta }$ vectors for the functions $q(t,1,3,\overline{\alpha },\overline{\beta },\gamma ),$ $q(t,2,3,\overline{\alpha },\overline{\beta },\gamma ),$ and $q(t,3,3,\overline{\alpha },\overline{\beta },\gamma )$ are shown in columns 1, 2, and 3, respectively.

$\mathit{q}(\mathit{t},1,3,\overline{\mathit{\alpha }}$,$\overline{\mathit{\beta }}$,$\mathit{\gamma }$)	$\mathit{q}(\mathit{t},2,3,\overline{\mathit{\alpha }}$,$\overline{\mathit{\beta }}$,$\mathit{\gamma }$)	$\mathit{q}(\mathit{t},3,3,\overline{\mathit{\alpha }}$,$\overline{\mathit{\beta }}$,$\mathit{\gamma }$)
${\alpha }_0={\beta }_0=-1.823·{10}^{-00}$	${\alpha }_0={\beta }_0=-2.384·{10}^{-00}$	${\alpha }_0={\beta }_0=8.617·{10}^{-01}$
${\alpha }_1={\beta }_1=2.598·{10}^{-03}$	${\alpha }_1={\beta }_1=7.752·{10}^{-03}$	${\alpha }_1={\beta }_1=6.852·{10}^{-03}$
${\beta }_2=-2.365·{10}^{-02}$	${\alpha }_2={\beta }_2=-1.191·{10}^{-06}$	${\alpha }_2={\beta }_2=-3.883·{10}^{-07}$
${\beta }_3=1.303·{10}^{-03}$	${\beta }_3=3.420·{10}^{-04}$	${\alpha }_3={\beta }_3=-4.681·{10}^{-11}$

for the functions $q(t,1,3,\overline{\alpha },\overline{\beta },\gamma ),$ $q(t,2,3,\overline{\alpha },\overline{\beta },\gamma ),$ and $q(t,3,3,\overline{\alpha },\overline{\beta },\gamma )$.

In addition, the family of candidate functions with $j=2$ was tested. The $RR$ value relating to the candidate function with $i=1$ and $j=2$ was high ($R{R}_{A2}\left(1\right)=0.9215$), while the residuals between experimental and predicted values oscillated between $-4.073$ and $7.253$ $\mathsf{\mu }$m at the experimental times (Table S10 and Figure S10). This function was able to show a good positive correlation between predictions and experimental data ($R^2=0.94$, Figure 5a), and the Shapiro–Francia normality test resulted in $p=0.84$ for in silico values (Figure 5b). Furthermore, an unpaired Student t-test between experimental data and in silico predictions resulted in $p=0.85$ (Figure 5c). Nevertheless, a plot of the function in the experimental range showed that no vertical asymptotes or stationary points were present (Figure 5d). Then, for this function $n_d=0$ and $n_z=0$.

The $RR$ value associated with the candidate function with $i=2$ and $j=2$ was elevated ($R{R}_{A2}\left(2\right)=0.9136$), while the residuals between experimental and predicted values oscillated between $-3.966$ and $11.714$ $\mathsf{\mu }$m (Table S11 and Figure S11). This function was able to show a good positive correlation between predictions and experimental data ($R^2=0.94$, Figure 5e), and the Shapiro–Francia normality test resulted in $p=0.16$ for in silico values (Figure 5f). Furthermore an unpaired Student t-test between experimental data and in silico predictions resulted in $p=0.86$ (Figure 5g). However, a complete plot revealed a stationary point at $t=4027.717$ h (Figure 5h). Therefore, for this function $n_d=0$ and $n_z=1$.

The $RR$ value referring to the candidate function with $i=3$ and $j=2$ was considerable ($R{R}_{A2}\left(3\right)=0.8982$), while the residuals between experimental and predicted values oscillated between $-1.593$ and $13.359$ $\mathsf{\mu }$m (Table S12 and Figure S12). This function was able to show a clear positive correlation between predictions and experimental data ($R^2=0.93$, Figure 5i), and the Shapiro–Francia normality test resulted in $p=0.11$ for in silico values (Figure 5j). Furthermore, an unpaired Student t-test between experimental data and in silico predictions resulted in $p=0.78$ (Figure 5k). The complete plot of the function showed a stationary point at $t=4112.742$ h (Figure 5l). Thus, for this function $n_d=0$ and $n_z=1$. The numeric values of the $\overline{\alpha }$ and $\overline{\beta }$ elements are detailed in

Table 4. The numeric values of the elements of $\overline{\alpha }$ and $\overline{\beta }$ vectors for the functions $q(t,1,2,\overline{\alpha },\overline{\beta },\gamma ),$ $q(t,2,2,\overline{\alpha },\overline{\beta },\gamma ),$ and $q(t,3,2,\overline{\alpha },\overline{\beta },\gamma )$ are shown in columns 1, 2, and 3, respectively.

$\mathit{q}(\mathit{t},1,2,\overline{\mathit{\alpha }}$,$\overline{\mathit{\beta }}$,$\mathit{\gamma }$)	$\mathit{q}(\mathit{t},2,2,\overline{\mathit{\alpha }}$,$\overline{\mathit{\beta }}$,$\mathit{\gamma }$)	$\mathit{q}(\mathit{t},3,2,\overline{\mathit{\alpha }}$,$\overline{\mathit{\beta }}$,$\mathit{\gamma }$)
${\alpha }_0={\beta }_0=7.707·{10}^{-01}$	${\alpha }_0={\beta }_0=8.033·{10}^{-01}$	${\alpha }_0={\beta }_0=8.469·{10}^{-01}$
${\alpha }_1={\beta }_1=2.449·{10}^{-03}$	${\alpha }_1={\beta }_1=1.86·{10}^{-02}$	${\alpha }_1={\beta }_1=1.058·{10}^{-02}$
${\beta }_2=2.469·{10}^{-03}$	${\alpha }_2={\beta }_2=-2.309·{10}^{-06}$	${\alpha }_2={\beta }_2=1.011·{10}^{-06}$
		${\alpha }_3=-4.032·{10}^{-10}$

for the functions $q(t,1,2,\overline{\alpha },\overline{\beta },\gamma ),$ $q(t,2,2,\overline{\alpha },\overline{\beta },\gamma ),$ and $q(t,3,2,\overline{\alpha },\overline{\beta },\gamma )$.

Moreover, the family of candidate functions with $j=1$ was tested. The $RR$ value related to the candidate function with $i=1$ and $j=1$ was considerable ($R{R}_{A1}\left(1\right)=0.88$), while the residuals between experimental and predicted values oscillated between $-11.781$ and $4.457$ $\mathsf{\mu }$m at the experimental times (Table S13 and Figure S13). This function was able to show a good positive correlation between predictions and experimental data ($R^2=0.88$, Figure 6a), and the Shapiro–Francia normality test resulted in $p=0.27$ for in silico values (Figure 6b). Furthermore an unpaired Student t-test between experimental data and in silico predictions resulted in $p=0.84$ (Figure 6c). In addition, a graphical representation of the function within the experimental range showed a stationary point at $t=3869.891$ h (Figure 6d), so, for this function, $n_d=0$ and $n_z=1$.

The $RR$ value associated with the candidate function with $i=2$ and $j=1$ was high ($R{R}_{A1}\left(2\right)=0.9093$), while the residuals between experimental and predicted values oscillated between $-3.728$ and $12.131$ $\mathsf{\mu }$m (Table S14 and Figure S14). This candidate function was able to show a good positive correlation between predictions and experimental data ($R^2=0.92$, Figure 6e), and the Shapiro–Francia normality test resulted in $p=0.15$ for in silico values (Figure 6f). Furthermore, an unpaired Student t-test between experimental data and in silico predictions resulted in $p=0.85$ (Figure 6g). However, a complete plot showed a stationary point at $t=4023.043$ h (Figure 6h). Therefore, for this function $n_d=0$ and $n_z=1$.

Similarly, the $RR$ value referring to the candidate function with $i=3$ and $j=1$ was considerable ($R{R}_{A1}\left(3\right)=0.899$), while the residuals between experimental and predicted values oscillated between $-1.307$ and $12.569$ $\mathsf{\mu }$m (Table S15 and Figure S15). This function was able to show a good positive correlation between predictions and experimental data ($R^2=0.90$, Figure 6i), and the Shapiro–Francia normality test resulted in $p=0.075$ for in silico values (Figure 6j). Furthermore, an unpaired Student t-test between experimental data and in silico predictions resulted in $p=0.85$ (Figure 6k). The complete plot of the function showed a stationary point at $t=4118.667$ h (Figure 6l). Thus, also for this function $n_d=0$ and $n_z=1$.

Finally, the $RR$ value pertaining to the candidate function with $i=4$ and $j=1$ was high ($R{R}_{A1}\left(4\right)=0.93$), while the residuals between experimental and predicted values oscillated between $-11.242$ and $3.514$ $\mathsf{\mu }$m (Table S16 and Figure S16). This function was able to show a good positive correlation between predictions and experimental data ($R^2=0.93$, Figure 6m), and the Shapiro–Francia normality test resulted in $p=0.44$ for in silico values (Figure 6n). Furthermore, an unpaired Student t-test between experimental data and in silico predictions resulted in $p=0.92$ (Figure 6o). However, the complete plot of the function showed three stationary points at $t=868.540$, $t=1677.503$, and $t=4445.434$ h (Figure 6p). As a consequence, for this function $n_d=0$ and $n_z=3$. The numeric values of the $\overline{\alpha }$ and $\overline{\beta }$ elements are listed in

Table 5. The numeric values of the elements of $\overline{\alpha }$ and $\overline{\beta }$ vectors for the functions $q(t,1,1,\overline{\alpha },\overline{\beta },\gamma ),$ $q(t,2,1,\overline{\alpha },\overline{\beta },\gamma ),$$q(t,3,1,\overline{\alpha },\overline{\beta },\gamma ),$ and $q(t,4,1,\overline{\alpha },\overline{\beta },\gamma )$ are shown in columns 1, 2, 3, and 4, respectively.

$\mathit{q}(\mathit{t},1,1,\overline{\mathit{\alpha }},\overline{\mathit{\beta }},\mathit{\gamma }$)	$\mathit{q}(\mathit{t},2,1,\overline{\mathit{\alpha }},\overline{\mathit{\beta }},\mathit{\gamma }$)	$\mathit{q}(\mathit{t},3,1,\overline{\mathit{\alpha }},\overline{\mathit{\beta }},\mathit{\gamma }$)	$\mathit{q}(\mathit{t},4,1,\overline{\mathit{\alpha }},\overline{\mathit{\beta }},\mathit{\gamma }$)
${\alpha }_0={\beta }_0=-1.047·{10}^{-03}$	${\alpha }_0={\beta }_0=7.879·{10}^{-01}$	${\alpha }_0={\beta }_0=8.575·{10}^{-01}$	${\alpha }_0={\beta }_0=4.432·{10}^{-01}$
${\alpha }_1={\beta }_1=-5.603·{10}^{-03}$	${\alpha }_1={\beta }_1=1.898·{10}^{-02}$	${\alpha }_1={\beta }_1=9.782·{10}^{-03}$	${\alpha }_1={\beta }_1=7.756·{10}^{-02}$
	${\alpha }_2=-2.46·{10}^{-06}$	${\alpha }_2=1.405·{10}^{-06}$	${\alpha }_2=-7.827·{10}^{-05}$
		${\alpha }_3=-4.457·{10}^{-10}$	${\alpha }_3=2.874·{10}^{-08}$
			${\alpha }_4=-3.09·{10}^{-12}$

for the functions $q(t,1,1,\overline{\alpha },\overline{\beta },\gamma ),$$q(t,2,1,\overline{\alpha },\overline{\beta },\gamma ),$$q(t,3,1,\overline{\alpha },\overline{\beta },\gamma ),$ and $q(t,4,1,\overline{\alpha },\overline{\beta },\gamma )$.

3.3. A Novel Metric to Discriminate Compatible Functional Forms

The previous analysis of each candidate function was able to link the standard statistics $R^2$ accounted for in the matrix $RR$ to the indexes $n_d$ and $n_z$. This knowledge was condensed in the novel matrix $\mathsf{\Xi }(RR,n_d,n_z)$ (Figure 7). This kind of matrix, which was formed by a 3$\times$5 matrix together with a singleton 1$\times$1 related to the combination $i=4$ and $j=1$, was able to account for the novel metric in Equation (16), which was used to assess the suitability of each candidate function. More specifically, this metric ranged between $-1.483$ for the candidate with $i=1$, $j=4$ and $0.9215$ for the candidate with $i=1$, $j=2$, while it was zero for the candidate functions with $i=1,2$ for the family with $j=3$. Similarly, it was zero both for $i=2,3$ and the families with $j=4$ and $j=5$. The families with positive values of the $\mathsf{\Xi }$ metric were $j=1$ and $j=2$, together with the combination $i=3$, $j=3$ and the singleton matrix (Figure 7).

3.4. Sensitivity Analysis of the Best Functional Form

The best functional form previously identified was $q(t,1,2,\overline{\alpha },\overline{\beta },\gamma )$. To further explore the influence of each parameter on the global evolution of this function, the increase in ${\alpha }_0$ through $\Delta {\alpha }_0$ was tested for changes in the curve shape (Figure 8a). In addition, the sensitivity of the candidate function to $\Delta {\alpha }_0$ (Equation (17)) was studied (Figure 8b) and resulted in a time-variant evolution increasing with $\Delta {\alpha }_0$ and decreasing with time. Similarly, the changes in shape for the increase in ${\alpha }_1$ through $\Delta {\alpha }_1$ were investigated (Figure 8c), while the sensitivity with respect to the increase in $\Delta {\alpha }_1$ (Equation (18)) is plotted in Figure 8d, and stabilised with time. Again, the changes in the curve shape for changes in ${\beta }_2$ through the increase in $\Delta {\beta }_2$ are explored in Figure 8e, and the sensitivity with respect to $\Delta {\beta }_2$ is plotted in Figure 8f and vanishes with the time increase. Finally, the shape changes of the best functional form through the increase in $\Delta \gamma$ are studied in Figure 8g, while the sensitivity with respect to $\Delta \gamma$ (Equation (19)) is plotted in Figure 8h and decreases with time.

3.5. Optimisation of Parameters of the Best Functional Form

Once the general influence of the parameters ${\alpha }_0,{\alpha }_1,{\beta }_2,\gamma$ on the best candidate function $q(t,1,2,\overline{\alpha },\overline{\beta },\gamma )$ was studied, their optimisation was performed. In particular, the ability of this function to model the experimental set of data was tested with respect to the variation in the standard $R^2$ statistic with the $\gamma$ parameter. The trend in the $R^2$ value decreased (from around $0.86$ to $0.75$) with the increase in $\gamma$ from around $5.65$ to $7.00$, except for the presence of sudden changes at around $5.82$ to $6.1$, and at around $6.5$ to $6.75$. In these intervals, the $R^2$ value increased beyond $0.95$ (see Figure 9a). More specifically, the maximum value ($R^2=0.9716$, $RMSE=3.247$) was reached for $\overline{\gamma }=6.5575$. For $\gamma =\overline{\gamma }$, the element of the parameter vectors was ${\alpha }_0=-1.201$, ${\alpha }_1=-2.004·{10}^{-3}$, and ${\beta }_2=-5.023·{10}^{-4}$. Finally, the modelling capabilities of the best candidate function $q(t,1,2,{\alpha }_0,{\alpha }_1,{\beta }_2,\overline{\gamma })$ were tested. This function resulted in a correlation between experimental data and predictions ($R^2=0.90$) for experimental determination times (Figure 9b), and it was also able to maintain the amount of residuals between $-14.388$ and $2.167$. The form of distribution of in silico predictions was compared to the form of distribution of experimental data and resulted in $p=0.3$ (Shapiro–Francia normality test $\alpha =0.05$, Figure 9c). The box plots of both experimental data and in silico predictions were compared (Figure 9d) through an unpaired Student t-test, which resulted in $p=0.77$. The time evolution of the best candidate function was, then, studied and resulted in a smooth and continuum function, which was plotted with $99.95\%$ confidence predictions bounds, and compared to the mean experimental data (mean values ± standard deviation). Finally, the increment in the $R^2$ standard statistics achieved in this work was calculated with reference to previous literature works [57,58] (Figure 9f).

4. Discussion

In this work, a computationally light and powerful approach was proposed to model the outgrowth of the scar tissue layer around implanted Parylene C devices over a time range of around 5400 h.

4.1. Some Restrictions Due to the Low Cardinality of the Experimental Set

The low cardinality of the experimental data set made the analysis of the potentially suitable candidate functions complex, because of the limited range of variation of both i and j indexes, which decreased the number of possible functional choices. More specifically, both indexes ranged between 1 and 5, since 1 is the minimum value for polynomials to provide a time-variant evolution, while 5 was previously known to be the minimum degree for polynomials to allow the onset of Runge instability, for a similar low-cardinality set ($n=9$) [57].

4.2. The Standard Optimisation Process

A standard optimisation was performed through a high-performance computational software (Matlab, R2022 64-bit, Academic, The MathWorks, Inc. USA) providing a selection between potentially suitable and not suitable combinations of the i and j indexes. The W and twin LRR matrices were used for this aim. Indeed, the W matrix provided the base 10 logarithm of the standard RMSE statistics, while the LRR matrix provided the base 10 logarithm of the standard $R^2$ statistics. Unsuitable candidates were recognisable through values greater than 1 of the W matrix elements (see Figure 1a). These values derived from a lack of convergence, because of the rise in Runge instability [61], which resulted in large numerical oscillations (e.g., all the functions with $i=5$ and the functions with $i=4$ and $2\le j\le 5$, see Figure 1b,c). The same nonlinear scale (base 10 logarithm) was applied to the twin LRR, which accounted for the standard statistics $R^2$. As a consequence, the same pattern of potentially suitable and not suitable combinations of indexes was transferred to the LRR matrix. Similarly, since this division between potentially suitable and not suitable candidate functions was invariant with respect to the change in scale, the same pattern was also applied to the $RR$ matrix, which was used as the main reference in this work.

4.3. The Standard Validation of Potentially Suitable Candidate Functions

The validation of potentially suitable candidates was performed by comparing their in silico predictions for optimised values of $\overline{\alpha }$, $\overline{\beta }$, and $\gamma$ (nonlinear least-squares method, least absolute residuals, Levemberg–Marquardt algorithm, see Table 1, Table 2, Table 3, Table 4 and Table 5) to the set of experimental data. In all cases, a quite strong positive correlation were found, since the standard $R^2$ statistics ranged between $0.88$ and $0.98$ for all potentially suitable combinations (see Figure 2a,e,i; Figure 3a,e,i; Figure 4a,e,i; Figure 5a,e,i; and Figure 6a,e,i,m). In addition, the statistical significance of the difference between in silico predictions and experimental data was tested. First, the form of the in silico distribution was explored through a Shapiro–Francia normality test ($\alpha =0.05$), which resulted in p values ranging from $0.07$ and $0.94$ (see Figure 2b,f,j; Figure 3b,f,j; Figure 4b,f,j; Figure 5b,f,j; Figure 6b,f,j,n) confirming the normal distribution of both in silico predictions and experimental data ($p=0.48$). Then, the significance of the difference between predictions and data was explored through an unpaired Student t-test, which resulted in p values ranging from $0.35$ to $0.99$ (see Figure 2c,g,k; Figure 3c,g,k; Figure 4c,g,k; Figure 5c,g,k; and Figure 6c,g,k,o), confirming the absence of any statistically significant difference between in silico predictions and experimental data. An explanation of this unexpected confounding phenomenon is related to the low cardinality of the experimental data set, due to the coarse time interval between experiments. In other words, several candidate functions, arising from different combinations, were able to replicate in a satisfactory way the experimental data just for selected experimental time points. Nevertheless, through this kind of analysis, no information was obtained about the behaviour of the candidate functions between two consecutive experimental time points.

4.4. Further Validation of Potentially Suitable Candidate Functions

Therefore, the exclusion of not suitable combinations of i and j through the standard procedure was not enough, because of the novelty of the proposed approach. As a consequence, each previously selected suitable candidate function was continuously plotted over the whole experimental time range. Through this further analysis, two different behaviours with inconsistent physical meaning (or with unclear physical meaning) were discovered: some candidate functions presented vertical asymptotes, while others had two or more stationary points along their evolution. The presence of vertical asymptotes was highly nonphysical, since it contradicted the continuity of the candidate function and predicted an infinite thickness of the scar tissue capsule. Similarly, the presence of two or more stationary points (i.e., points where the first derivative of the thickness over time with respect to the time was zero) was likely due to numeric oscillations, since alternate phases of thickness increase and decrease were not supported by experimental data.

4.5. The Novel $\mathsf{\Xi }$ Metric

The presence of nonphysical behaviours was accounted for in the novel metric expressed through Equation (16), where the numbers of vertical asymptotes and of stationary points were both considered for each candidate function to offset the apparent high values of the $RR$ matrix. Through this novel metric, a different picture of the global suitability of each candidate function arose (see Figure 7). More specifically, this matrix was represented by a 3$\times$5 matrix together with a singleton matrix for the combination $i=4$, $j=1$. In particular, on the right side of the matrix were grouped some combinations, which were unsuitable because of the presence of one ore more vertical asymptotes in the continuum evolution of the candidate function. In particular, the candidate functions resulting from the combinations $i=1$, $j=4$, and $i=1$, $j=5$ presented three and two vertical asymptotes, respectively. Differently, the candidate with $i=3$ and $j=3$ presented a low value of the RR matrix, while the combination $i=4$, $j=1$ resulted in numerical instability due to the Runge phenomenon. On the contrary, the left side of the $\mathsf{\Xi }$ matrix grouped some potentially suitable combinations, which were all compatible with the physical evolution of the thickness outgrowth, and among them was the best candidate function, for the combination $i=1,j=2$.

4.6. Sensitivity and Optimisation of the Best Candidate Function

Both the shape change and the sensitivity of the best functional form $q(t,1,2,\overline{\alpha },\overline{\beta },\gamma )$ were then explored with respect to numerical changes in relevant parameters belonging to the $\overline{\alpha }$ and $\overline{\beta }$ vectors (i.e., ${\alpha }_0,{\alpha }_1,{\beta }_2$), as well as with respect to the $\gamma$ value (see Figure 8). First, the change in shape of $q(t,1,2,\overline{\alpha },\overline{\beta },\gamma )$ for increasing values of $\Delta {\alpha }_0$ was studied and resulted in a shift in the curve towards higher values, while the sensitivity of the candidate function started from 1 and decreased over time. Similarly, for increments in $\Delta {\alpha }_1$ the shape of the curve was progressively more steep, while the sensitivity with respect to these changes stabilised over time towards an increasing steady value near to 1. On the contrary, the shape of $q(t,1,2,\overline{\alpha },\overline{\beta },\gamma )$ was practically unaffected by increasing changes in $\Delta {\beta }_2$, while the sensitivity with respect to theses changes started from around $0.1$ and decreased over time. Finally, the changes in $\Delta \gamma$ affected the curvature of the best candidate function, as well as the sensitivity of the curve, which started from values increasing with $\Delta \gamma$ and decreased over time. As a consequence, this analysis was able to show that the influence of the ${\alpha }_1$ parameter dominated the shape changes of the best candidate function, while the influence of $\gamma$ was greater than the influence of ${\alpha }_0$, which is, in its turn, greater than the influence ${\beta }_2$. In addition, the value of the standard statistics $R^2$ was measured for different values of the $\gamma$ parameter. A global decreasing trend was found together with a second, more complex, transition phenomenon. Indeed, the value of $R^2$ suddenly increased by $0.144$ when $\gamma$ increased by $0.001$, while it suddenly decreased by around $0.13$ when $\gamma$ increased by $0.05$. Similarly, a second similar transition was found, where the $R^2$ value suddenly increased by $0.1$ when $\gamma$ increased by $0.001$ and decreased by $0.2$ for an increment of $0.3$. Moreover, the value of $R^2$ was almost constant (even if slightly decreasing) for two numeric intervals (see Figure 9a). This novel kind of numeric instability was likely related to a boundary condition problem. Indeed, the candidate functions should not only model the experimental behaviour along the experimental time range, but they should also respect the initial boundary conditions. The proposed framework resulted in nonlinear initial boundary conditions; thus, standard numeric optimisation procedures, which provide algebraical approximations of transcendent functions, could become unstable around given values of the $\gamma$ parameters, creating a sudden numerical transition of the $R^2$ values. Therefore, the values of ${\alpha }_0,{\alpha }_1,{\beta }_2$ and $\gamma$ were optimised (nonlinear least-squares method, least absolute residuals, Levemberg–Marquardt algorithm) in order to keep the best candidate function within the stability range of the $\gamma$ parameter. This best candidate function was able to provide a good positive correlation with experimental data (see Figure 9b). In addition, the form of the distribution of in silico predictions was normal ($p=0.3$, Shapiro–Francia normality test with $\alpha =0.05$), and it was not statistically different from experimental data ($p=0.77$ unpaired Student t-test). The magnitude of the p value of this test was so high that any difference was excluded for every standard level of statistical significance (e.g., $\alpha =0.05$, as well as $\alpha =0.01$). Moreover, the continuum evolution of the best candidate function was able to provide a good modelling of the mean experimental data $R^2=0.97$, while the $99.95\%$ confidence prediction bounds were able to include the data error bars (see Figure 9e).

5. Conclusions

In this work, a novel theoretical framework was proposed to model the evolution of the capsule thickness over time together with an improved optimisation procedure able to avoid apparently suitable choices, as provided by standard numeric procedures. The proposed framework resulted in an fast, simple, and accurate modelling of experimental data ($R^2=0.97$), definitely improving on previous literature approaches [57,58] ($R^2=0.86$ and $R^2=0.92$, respectively (see Figure 9f)). Although the proposed approach was experimentally validated for rat peripheral nerves and inserted Parylene C devices, it may be, more generally, able to model the evolution of the scar tissue thickness around implanted polymeric structures in vivo. Indeed, the formulation of this novel phenomenological approach was not related to the specificity of the tissue, while the selected suitable (compatible and best) functional forms were able to considerably change their shape for changes in values of vectors $\overline{\alpha },\overline{\beta }$. Moreover, the selected functional forms were also able to be superimposed to model more complex behaviours. As a consequence, it could be used as a benchmark to compare the immune response of peripheral nerves to more complex devices, involving both mechanical and electrical stimuli and resulting in scar tissue production. Similarly, once the relationship between scar tissue thickness and the impedance of devices is known, it may also be used to assess their possible lifetime, in order to prevent their electrical failure.