bioengineering Article Viscoelastic Characterization of Parasagittal Bridging Veins and Implications for Traumatic Brain Injury: A Pilot Study Silvia García-Vilana 1,† , David Sánchez-Molina 1, *,† , Jordi Llumà 1 , Ignasi Galtés 2,3 , Juan Velázquez-Ameijide 1 , M. Carmen Rebollo-Soria 4 and Carlos Arregui-Dalmases 5 1 Escola d’Enginyeria de Barcelona Est, Universitat Politècnica de Catalunya, Av. Eduard Maristany, 16, 08019 Barcelona, Spain; silvia.garcia.vilana@upc.edu (S.G.-V.); jordi.lluma@upc.edu (J.L.); juan.velazquez@upc.edu (J.V.-A.) 2 Institut de Medicina Legal i Ciències Forenses de Catalunya, G.V. Corts Catalanes, 111, 08014 Barcelona, Spain; ignasigaltes@gmail.com 3 Departament de Psiquiatria i de Medicina Legal, Universitat Autònoma de Barcelona, Cerdanyola del Vallès, 08193 Barcelona, Spain 4 Instituto de Medicina Legal de Aragón, Irene Izárbez, 22005 Huesca, Spain; mcrebollos@aragon.es 5 Centro Zaragoza, Crtra. 232, km.273, 50690 Pedrola, Spain; carlos.arregui@centro-zaragoza.com * Correspondence: david.sanchez-molina@upc.edu † These authors contributed equally to this work. Abstract: Many previous studies on the mechanical properties of Parasagittal Bridging Veins (PSBVs) found that strain rate had a significant effect on some mechanical properties, but did not extensively   study the viscoelastic effects, which are difficult to detect with uniaxial simple tensile tests. In this study, relaxation tests and tests under cyclic loading were performed, and it was found that PSBVs Citation: García-Vilana, S.; do indeed exhibit clear viscoelastic effects. In addition, a complete viscoelastic model for the PSBVs is Sánchez-Molina, D.; Llumà, J.; Galtés, I.; proposed and data from relaxation, cyclic load and load-unload tests for triangular loads are used to Velázquez-Ameijide, J.; Rebollo-Soria, M.C.; Arregui-Dalmases, C. find reference values that characterize the viscoelastic behavior of the PSBVs. Although such models Viscoelastic Characterization of have been proposed for other types of blood vessels, this is the first study that clearly demonstrates Parasagittal Bridging Veins and the existence of viscoelastic effects from an experimental point of view and also proposes a specific Implications for Traumatic Brain model to explain the data obtained. Finally, this study provides reference values for the usual Injury: A Pilot Study. Bioengineering viscoelastic properties, which would allow more accurate numerical simulation of PSBVs by means 2021, 8, 145. https://doi.org/ of computational models. 10.3390/bioengineering8100145 Keywords: bridging veins; TBI; tissue characterization; biomechanics; strain rate dependent materi- Academic Editor: Liang Luo als; viscoelasticity Received: 19 August 2021 Accepted: 12 October 2021 Published: 18 October 2021 1. Introduction Publisher’s Note: MDPI stays neutral The importance of the TBI at the global level cannot be overlooked. Every year, about with regard to jurisdictional claims in 70 million people worldwide are afflicted with some form of TBI, with the global incidence published maps and institutional affil- reaching about 940 cases per 100,000 people [1]. Among the factors leading to different iations. types of TBI, a particularly pressing one is SDH produced by the rupture of some types of PSBVs [2,3]. In particular, the superior sagittal sinus-bridging veins, which link the central sagittal sinus to the upper part of the brain mass, have been identified as a highly critical region prone to producing SDHs of some severity [4–7]. Copyright: © 2021 by the authors. For this reason, detailed and accurate studies of the mechanical properties of the Licensee MDPI, Basel, Switzerland. collagenous tissue constituting the PSBVs are important, and should be updated using This article is an open access article better and more accurate measurement techniques. For example, one influential paper distributed under the terms and detected a negative correlation of ultimate stress and strain rate [4]; however, a later paper conditions of the Creative Commons showed a positive correlation [7], although neither author found the correlation to be Attribution (CC BY) license (https:// statistically significant. The issue was recently closed when the use of more accurate creativecommons.org/licenses/by/ digital measurement systems found a positive and significant correlation [3]. Although 4.0/). strain-rate dependent behavior for some mechanical properties have been studied by Bioengineering 2021, 8, 145. https://doi.org/10.3390/bioengineering8100145 https://www.mdpi.com/journal/bioengineering Bioengineering 2021, 8, 145 2 of 15 several authors [4,6–8], there are still few studies in the literature properly investigating the viscoelastic behavior of PSBVs, including the measurement of relaxation times or preconditioning under cyclic loading [3,9]. An important advance is found in [10] which examines preconditioning by cyclic loading and fatigue. However, although some work has been done for other types of blood vessels [11], no specific and complete viscoelastic model has yet been proposed for the PSBVs. The aim of this study is to contribute to filling the gap in viscoelastic characterization, using a more complex series of tests including load-unload cycles, stress relaxation cycles and fast repeated loading cycles, following the line of other works on the characterization of viscoelastic behavior of biological tissues [12–14]. This work shows that, the viscoelastic behavior of the PSBVs is clearly observable, measurable and quantifiable when observing the rehology of the tissue in a series of more complex tests, such as cyclic-load and relaxation tests, instead of conventional simple uniaxial tests, where the contribution of the viscoelastic effect to the mechanical behavior is difficult to quantify [3]. Such viscoelastic behavior can be used to improve both the understanding of injury mechanisms [15,16], and the computational models used to assess the probability of severe injuries. 2. Data and Methods 2.1. Material and Specimen Preparation For this study, different sections of the meningeal-cortex space were obtained from four autopsies of post-mortem human subjects (PMHS) conducted at the Forensic Pathology Service of the Legal Medicine and Forensic Science Institute of Catalonia (FPS/IMLCFC). The study was approved by the Research Committee of the IMLCFC. None of the subjects had been previously diagnosed with any blood vessel pathology. Once received, the meningeal sections were kept refrigerated for, at most, 96 h in airtight containers that maintained the natural degree of hydration. The sections of the meningeal-cortex space included the meninges, the subarachnoid space, and the upper part of the cerebral cortex. From those four sections, six PSBVs were carefully dissected (more than one PSBV was dissected from two sections). After dissection, the PSBVs were kept in contact with towel tissue soaked in phosphate-buffered saline solution. Twelve to twenty-four hours before testing, the PSBVs were placed in a refrigerator at 2 ◦ C and before testing they were allowed to come to room temperature for one hour. 2.2. Cyclic and Relaxation Tests The tests were performed with an Allround-Table-Top Universal Test Machine (UTM) Zwick-Roell® and load was measured with a 20N load cell HBM® . Special fixtures were used for the PSBVs tests and the displacement was carefully measured with a digital control unit attached to the UTM (due to the small dimensions of the cross-section of the PSBVs, no significant local effects appear in the displacement measures). The specimens were kept in physiological saline solution at 5 ◦ C to minimize tissue degradation during the testing. During the tests, a digitally displacement-controlled maximum strain ε max was applied and the force experienced by each PSBV specimen was measured very accurately. In order not to cause permanent damage to the specimens, ε max was selected as 20% of the average ultimate strain, determined in previous tensile studies of PSBVs [3]. The sequence of stages of the loading process was similar to that used in other soft tissue studies such as [12]. The different stages are as follows: (1) load-hold at ε max , (2) rapid load-unload cycles at ε max , (3) three stages of load-hold at ε max , 2/3·ε max and 1/3·ε max , (4) three load-unload stages at ε max at stretch rates of 0.01, 0.1 and 1 s−1 and (5) load to failure. Figure 1b shows a schematic of the above stages. The tests were performed at a speed of 50 mm/s (except for stage 4) and the holding time and the time between load stages were 50 s, allowing the specimen to recover completely. Bioengineering 2021, 8, 145 3 of 15 (a) (b) Figure 1. (a) Experimental setting measurements: the coordinate system is centered on the static fixture with the x axis parallel to the load direction, and a position change from the initial X to the deformed xt configuration is represented. (b) Cyclic and relaxation test stages scheme: (1) load-holding to maximum strain (ε max ), (2) fast load/unload cycle stage to maximum strain, (3) load-holding to ε max , 2/3·ε max and 1/3·ε max , (4) load and unload to ε max at 0.01, 0.1, 1 s−1 and (5) load to fracture. In this study, each PSBV is modeled as a homogeneous, transversely isotropic, hollow cylinder with constant cross-section [3–5,7]. The coordinate axes have their origin upon the static bottom fixture; the X axis is aligned with the (longitudinal) force, while the Y and Z axes are contained in the PSBV cross-section. Thus, the stretch for each time instant is easily determined as λt = `t /`0 = 1 + δt /`0 where `0 is the initial length, `t the instantaneous length and the displacement δt = `t − `0 has been digitally measured. Let X t = ( Xt , Yt , Zt ) be the coordinates in the initial configuration and xt = ( xt , yt , zt ) in the deformed configuration, the deformation of the cylinder is expressed as h i1/2 h i1/2 xt = Xλt , yt = Y 1 − ν̄(λ2t − 1) , zt = Z 1 − ν̄(λ2t − 1) (1) where ν̄(λt ) represents the Poisson effect. Thus, being Ft = ∂xt /∂X the deformation gradient tensor, the (Green-Lagrangian) strain tensor ε = (FtT Ft − 1)/2 is obtained as λ2 − 1   0 0 1 t ε(t) = 0 −ν̄(λt )(λ2t − 1) 0  (2) 2 2 0 0 −ν̄(λt )(λt − 1) On the other hand, the force Ft has been measured for each time instant during the whole test. From the force, the (second Piola-Kirchhoff) stress tensor σ is obtained:   Ft λ A 0 0 σ (t) =  t 0  (3)  0 0 0 0 0 0 being A0 the initial (undeformed) cross-section of the PSBV. 2.3. Quasi-Linear Viscoelastic Model Many previous studies had used uniaxial tensile tests on PSBVs, showing that an elastic model reasonably explains the results obtained in this type of test. For this reason, the collagenous tissue of PSBVs has been treated simply as an elastic material. However, Bioengineering 2021, 8, 145 4 of 15 the cycle tests and the relaxation tests performed show that there is a clear viscoelastic effect on the mechanical behavior of PSBVs and that, therefore, viscoelastic models would work better for more general situations than simple traction. In this study, a viscoelastic model of type QLVE [13,17] will be used to model the data. In the QLVE model, the mechanical response is divided into a strain-rate dependent (e) (v) part and a strain-rate independent part: σx = σx + σx , using a separable relaxation (e) function R(ε x , t) = G (t) · ∂σx /∂ε x [13,18], where G (t) is the relaxation function, given by the well-known Prony series [19–22]: N G (t) = 1 + ∑ gk e−t/τk (4) k =1 with gk being the “weights” of each k-term and τk being the relaxation times, both to be obtained from experimental trials. The value of N depends on the needs of the setting, in Section 3.1 it is shown that N = 3 is sufficient to represent the data in this study. With these assumptions, the axial strain (using the second Piola-Kirchhoff stress tensor) is expressed as (e) ∂Ψ Z t ∂σ σx (t) = = G (t − τ ) x ε̇ x (τ ) dτ (5) ∂ε x 0 ∂ε x 3. Results For this study, cyclic and uniaxial relaxation tests were performed on six PSBV speci- mens. Figure 2 clarifies the five stages of the testing process for each specimen. The same figure shows the variation of the imposed displacements over time, together with the axial force measured by the load cell. Figure 2. Imposed displacements in the different test stages and reaction force of the PSBV for specimen 2632A/20. The stages with constant displacement have been enlarged to show the decrease in force in stages (1) to (3), and the concavity of the load–discharge curve of stage (4). Stages (1) and (3) are stages with constant strain: the first one up to the established strain level ε max , and the second one divided in three stages up to ε max , 2/3 · ε max and Bioengineering 2021, 8, 145 5 of 15 1/3 · ε max . In both stages, while keeping constant the applied displacement and, therefore, the strain, a clear relaxation of the reaction force of the PSBV specimen is observed, with a drop of exponential type, which is a clear indication of the presence of viscoelastic effects in the mechanical behavior of the specimen. Likewise, stage (2) in which successive loading and unloading cycles have been applied up to ε max at high speed, the same behavior is perceived, showing a progressive drop in the force as the number of accumulated cycles increases. Finally, stage (4) is divided into three load-unload stages at speeds of 0.01, 0.1 and 1 s−1 respectively. It can be seen that, the loads being triangular and applied at a constant velocity, the force response is clearly a concave curve, which also indicates the presence of viscoelasticity. Thus, the results obtained from the cycling and relaxation tests indicate that the bridging vein is a viscoelastic material; therefore, a constitutive model that contemplates such viscoelasticity is required to describe its mechanical behavior. 3.1. Relaxation Tests The parameters of the viscoelastic model of Equation (5) were fitted from the data of stages (1) and (3) of Figure 2. For this purpose, the set of force-time values from the instant at which the maximum displacement is reached until the end of the stage was considered. The axial stress (σx ) is expressed in terms of the force as σx (t) = Ft /λt A0 , see Equa- tion (3). Given the low level of deformation, (ε < 5%), the material can be considered linear elastic [23] without making large errors and employing the following relationship: " # N Z t Ft = λ t ε x ( t ) + ∑ gk e−(t−τ )/τk ε̇ x (τ ) dτ (6) EA0 k =1 0 Note that, during the initial instants of loading, ε̇ x 6= 0 and that makes the integral not identically zero (see that once the maximum displacement ε̇ x = 0 is reached, the integral term is null). However, at the end of the stage, the specimen is completely relaxed, so the viscoelastic contribution (VC) ends up being zero. Moreover, from practically the beginning of the trial, both λt and ε x are constants: λt = λ∞ and ε x (t) = ε max . All this implies that the quotient Ft /EA0 also tends to the constant given by λ∞ ε max . Thus, measuring the force F∞ at the end of the stage provides a way to determine the value of EA0 from the ratio EA0 = F∞ /(λ∞ ε max ). For that reason, calculations can be made without a direct measurement of the elastic modulus or cross-section of the specimen by simply using the normalization f t = Ft /EA0 , where the latter magnitude f t will be called reduced force. During the relaxation itself, λt and, therefore, also ε x are constant. However, be- fore reaching their constant value, there is a very fast sudden stretching, during which the elongation is given by λt = 1 + δt /`0 , with δt being the displacement and the stretching rate being λ̇t = δ̇t /`0 . Experimentally, it was observed that the fast stretching velocity process can be adequately described by a function of “gamma distribution” type [24], with the form !n ∆max t e−t/ta δ̇t = (7) ta t a Γ ( n + 1) where t a is a parameter related to the time required to reach the maximum displacement and, in practice, is obtained by fitting the measured displacement velocity data. Similarly, several preliminary tests showed that taking n = 3 is sufficient to adequately represent the velocity at the beginning of each relaxation stage. Therefore, the displacement δt can be obtained by integrating the above expression (see Figure 3). Bioengineering 2021, 8, 145 6 of 15 Figure 3. Experimental displacement and experimental velocity applied to a specimen, plotted together with the fit by gamma-type function. The plot also shows the derivatives in Sections 1 and 3. It can be seen that the fit is in both cases more extreme than 0.997. From the explicit form (7) of the velocity δ̇t , we obtain the strain rate as ε̇ x (t) = δt which, introduced in Equation (6), leads to the formula: " # Ft N e−α1 t/θk 11 = λt ε max + ∑ gk (1 − θk )4 (1 − 2θk )7 m∑ ft = θkm (µm + νm e−α2 /θk ) (8) EA0 k =1 =1 In this formula, α1 and α2 are real numbers that depend on the initial conditions of the specimen (length) and the test velocity curve (time of load application, etc.). In the same way, the explicit analytical integration yields the real numbers µm and νm (the same calculation shows that νm > 0, µm < 0 for even m and νm < 0, µm > 0 if m is odd, and furthermore one has, µ1 , µ2 , µ3 , µ4 = 0). The only adjustable parameters in the above expression are θk = τk /t a and gk which are dimensionless. The relaxation times are obtained as τk = θk t a . The number of adjustable parameters is 2N, where N is the number of terms in the Prony series, see Equation (4). Figure 4 shows the measured force data for stage 1 of one of the specimens, and the corresponding fits with the viscoelastic model for N = 1, N = 2 and N = 3 in the Prony series (4). The fitted parameters have been set out in Table 1. It can be seen that as N increases, the fit improves markedly as expected, since the number of adjustable gk and gk parameters increases: 1. Notice that, although for N = 1 the coefficient r2 = 0.889 is quantitatively adequate, qualitatively it is observed that the fit does not adequately represent the vertical asymptote of force drop, nor the curve in general. 2. With N = 2 there is a marked improvement in the fit (r2 = 0.990). 3. Finally, for N = 3 (r2 = 0.995), both the asymptote and the curve are adequately represented by the viscoelastic model. Bioengineering 2021, 8, 145 7 of 15 Figure 4. Fitting of the proposed viscoelastic model, for stage (1) of one of the specimens. Fits with N = 1, N = 2 and N = 3 terms in the Prony series are considered, with the corresponding parameters from Table 1. Table 1. Parameters gk , τk of the viscoelastic model obtained for the specimen 635/21B for different values of k. N g1 [–] g2 [–] g3 [–] τ1 [s] τ2 [s] τ3 [s] N=1 0.360 — — 10.756 — — N=2 0.243 0.393 — 13.397 0.799 — N=3 0.215 0.221 0.398 18.162 1.885 0.211 Thus, in this study the parameters gk and τk have been determined for N = 3 ( g1 , g2 , g3 ; τ1 , τ2 , τ3 ). The viscoelastic model fits have been repeated for stages (1) and (3), the latter divided into subplots (3.1), (3.2) and (3.3), thus obtaining the parameters. The average for each specimen is presented in Table 2. Bioengineering 2021, 8, 145 8 of 15 Table 2. Average ± standard deviation of the parameters g1 , g2 , g3 and τ1 , τ2 , τ3 of the viscoelastic model in the four relaxation stages obtained for each specimen. Specimen g1 [–] g2 [–] g3 [–] τ1 [s] τ2 [s] τ3 [s] 2632A 0.10 ± 0.05 0.13 ± 0.04 1.70 ± 0.29 10.96 ± 2.64 0.88 ± 0.19 0.07 ± 0.01 617A 0.21 ± 0.10 0.32 ± 0.16 0.34 ± 0.06 15.38 ± 4.90 1.99 ± 2.62 0.17 ± 0.19 617B 0.24 ± 0.07 0.13 ± 0.08 0.30 ± 0.23 13.78 ± 9.15 0.96 ± 0.70 0.09 ± 0.08 621A 0.20 ± 0.07 0.18 ± 0.10 0.48 ± 0.42 17.08 ± 16.11 1.20 ± 1.23 0.09 ± 0.10 635A 0.12 ± 0.07 0.13 ± 0.07 0.84 ± 1.13 20.73 ± 13.34 2.16 ± 1.46 0.20 ± 0.10 635B 0.14 ± 0.06 0.16 ± 0.08 1.60 ± 1.62 25.41 ± 8.22 1.76 ± 1.93 0.17 ± 0.20 It can be observed in the table that the higher the value of gk (from g1 to g3 ), the lower the value of τk , thus increasing the degree of contribution of the model term, and the shorter the relaxation time. The relaxation times decrease by an order of magnitude progressively between τ1 and τ3 , obtaining similar values among the different specimens. The aver- age characteristic times for the sample are τ1 = 17.22 ± 5.17 s, τ2 = 1.49 ± 0.55 s and τ3 = 0.13 ± 0.05 s. The average values of the contribution coefficients for the sample are g1 = 0.17 ± 0.06, g2 = 0.17 ± 0.08 and g3 = 0.87 ± 0.62. In addition, to show the importance of the viscoelastic effect in PSBVs, the VC for the reduced force has been calculated, i.e., what percentage of the response corresponds to the elastic effect and what amount to the viscoelastic one. The VC for all specimens, calculated from the parameters gk and τk in Table 2, is shown in Figure 5. It can be seen that at the initial instant, the VC is between 26 and 35%, decreasing exponentially until it reaches approximately 0% at the end of the relaxation stretch. In fact, after 15 s of relaxation, the specimen has relaxed to a great extent and the VC is less than 10% in all samples. Figure 5. Evolution of VC in strength over time for the analyzed specimens: the black lines indicate the average elastic and viscoelastic contributions as a function of time, the shaded area indicates the error interval in the sample. Likewise, the average VC and viscoelasticity and the confidence interval have been determined from the parameters of the studied specimens (Figure 5). Thus, it can be seen that initially the viscoelastic part corresponds to 29.8 ± 3.35%, decreasing to 11.4 ± 2.40% at 10 s and 4.8 ± 1.06% at 25 s, while the elastic contribution starts from 69.8 ± 4.03%, rising to 88.5 ± 2.64% at 10 s to 94.9 ± 1.10% at 25 s. 3.2. Fast Loading/Unloading Cycle Tests The viscoelastic effect can also be analyzed in test stage (2), in which fast loading and unloading cycles have been performed. These cycles are fast compared to some of the relaxation times τk of the material. This leads to the assumption that in the time in Bioengineering 2021, 8, 145 9 of 15 which a loading and unloading cycle elapses, the BV may not have enough time to relax the internal stresses and this is reflected in the maximum forces reached in each cycle. Figure 6 shows the force peaks in successive cycles. It can be seen how in each new cycle the maximum force reached is lower, showing a progressive relaxation of the tissue of the PSBV as the number of cycles increases. For each sample, the USF (unconditioned scale factor), defined in [12], has been calculated as the following ratio: F1 USF(n, T ) = >1 (9) Fn where n is the cumulative cycle number and T is the period of the cycles (T = 0.55). As can be seen in the figure above, the magnitude of the USF grows progressively with the number of cycles. In some samples, it has reached 20 cycles, while in others it has been limited to 10 cycles. From the results, it has been calculated that the USF(n = 10) = 1.14 ± 0.085 and USF(n = 20) = 1.23 ± 0.095. In Appendix A, it is justified that the USF should be able to be approximated by a function of the type 1 + g0 USF(n, T ) ≈ (10) 1 + g0 e − α ( T ) n Figure 6. Force peaks in stage (2) as a function of the number of cycles performed and USF computed and predicted by Equation (10). This equation has the form of Equation (A7) and that prediction is just what is observed in Figure 6, where the fit between the above formula and the experimental data is very good (r2 = 0.983). The results obtained for the parameters in the sample set are g0 = 0.206 ± 0.047 and α = 0.127 ± 0.028. 3.3. Load-Unload Tests In addition to relaxation tests and fast loading cycle tests, three additional loading and unloading tests were performed in the form of triangular waves at very different strain rates, corresponding to stages 4a, 4b and 4c performed at constant stretching speeds λ̇t = 0.01, 0.1 and 1 s−1 . The main purpose was to verify whether the characteristic times (τk ) and contribution coefficients (gk ) obtained from the relaxation tests could well explain the shape of these curves. The displacement had a very roughly triangular shape (except for brief periods of acceleration and deceleration), each triangular wave included about 350 measurements along the loading and unloading cycle (Np = 285, 365 and 385 points, Bioengineering 2021, 8, 145 10 of 15 respectively), that allowed making a numerical integration of the viscoelastic contribution in (6), using Simpson’s 3/8 rule [25]: f (τ ) Z t = Np h z }| { Ik := eτ/τk ε̇ x (τ ) dτ ≈ 0 (11) N p −3 1 3h = 4 ∑ 8 ( f (tk ) + 3 f (tk+1 ) + 3 f (tk+2 ) + f (tk+3 )) k =0 where h = t/Np and tk = kh. The result can be seen in Figure 7 where the measured and model-predicted force curves are shown. The fit is reasonable although the medium and low velocity curves show larger divergences from the model. This could be partly due to the fact that the integration error is larger for the slower loading curves. This happens because the number of points is similar in the three curves, but the duration Tw of the triangular wave is different, the numerical error depends on the time step h used for the numerical integration, which is given by the relation h = Tw /Np (Np being similar for the three curves). The error resulting from Simpson’s 3/8 rule of integration is given by: 3h5 d4 f (ξ k ) 3Tw5 e≤ Np ∑ max < C ε̇ (12) 80 k ξ k ∈(tk−1 ,tk ) dτ 4 80Np4 Figure 7. Experimental force–time curves compared with model predictions: (a) fast curve (ε̇ ≈ 1 s−1 ), (b) medium velocity curve (ε̇ ≈ 0.1 s−1 ) and (c) slow curve (ε̇ ≈ 0.01 s−1 ). The predicted curve has been obtained by numerical integration. The cumulative integration error is larger for the slow curve and smaller for the fast curve. However, the fast curve seems to have a larger measurement error. Possibly, for that reason, the best fit is obtained for the medium speed curve. So, for the faster curve a smaller error was made, while for the slower curve the integration error was larger (having a similar number of points, but the duration of the last wave being longer, which makes the division coarser h = T/Np ). 4. Discussion The mechanical behavior of PSBVs has been previously studied by some researchers [5,7,9,26], although most of the previous studies do not explicitly consider viscoelastic effects in the stress/force–strain curves [10]. Some strain rate related effects are reported, especially with respect to the influence of strain rate on the mechanical failure of PSBVs [3,4,6,7]. Interestingly, like other biological tissues, PSBVs appear to exhibit failure stresses negatively correlated with strain rate [10,27]. However, the relaxation and fast-cycle tests performed in this study clearly show a viscoelastic effect in the PSBV response, especially in stress relaxation stages under constant strain, where the force on the PSBVs is progressively reduced. This fact had not been investigated in detail in the literature on mechanical properties of PSBV [9], although viscoelastic effects had been measured in rapid load cycles [10]. Bioengineering 2021, 8, 145 11 of 15 Regarding the results obtained, this work is possibly the first study of PSBVs estimat- ing the necessary relaxation times of the Prony series, see Equation (4), for modeling the vis- coelastic behavior of PSBVs. Interestingly, excellent fits (r2 > 0.99) are obtained with three relaxation times and these differ by an order of magnitude, as τ1 = O(101 ), τ2 = O(100 ), τ3 = O(10−1 ), so each of them seems to capture a different time scale. This is in agreement with the findings of Funk et al. (2000) [12], where the viscoelasticity of ankle ligaments is studied and where the observed scales are τ1 = O(102 ), τ2 = O(101 ), τ3 = O(100 ). On the other hand, in Davis & De Vita (2012) [13], it is pointed out that the QLVE model applied to rat tendons only captures the behavior well for t < 10 s, which is why they propose a nonlinear viscoelastic constitutive model that overcomes the shortcomings of QLVE to explain their data. As for the analysis of the preconditioning, which is based on the USF given by Equation (9), an analysis based on viscoelastic behavior has been presented here, which allows predicting that this magnitude will vary from one cycle to another by means of Equation (10). With only two parameters, this formula predicts a better fit than the exponential type heuristic formula: 3 USFn = a0 + ∑ a k e − bk n (13) k =1 used in [10] and previously in [28–30]. Although the Formula (13) is numerically adequate, it is a heuristic formula and is not directly derived from the viscoelastic behavior equations unlike the Formula (10). As for the applications of this work, they can be grouped into three different areas: 1. Clinicians are charged with the significant task of distinguishing between accidental and inflicted head trauma. Some times this distinction is straightforward, but in many cases the probabilities of injuries from accidental scenarios are unknown, making the differential diagnosis difficult [31]. A refinement of the knowledge of the tolerance ranges against rupture of PSBVs may provide greater accuracy in the reconstruction of injury mechanisms. 2. Computational biomechanics can simulate many potentially traumatic situations, so that models already allow detailed reconstructions of the sequence of events leading to a severe SDH. Accurate knowledge of the material behavior can improve FEHMs, as their inaccuracy is often not so much a computational problem, but a poor char- acterization of the biomechanical behavior of brain structures. For example, a good number of FEHMs use a stress-strain response for PSBVs that does not reflect the measured nonlinear behavior [3], e.g., the UDS FEHM (Université de Strasbourg) [32], the KTH FEHM (KTH Royal Institute of Technology) [33], UCDBTM (University College Dublin) [34,35], WSUBIM (Wayne State University) [36] or G/LHM [37] also model PSBVs as elastic beams with a linear stress-strain response [9]. The recognition of the importance of the nonlinear behavior and viscoelasticity of brain structures has been explicitly pointed out by the developers of YEAHM (University of Aveiro) [38] and interestingly, some FEHMs model PSBVs as nonlinear elastic materials [39]. In particular, Equation (4), together with the averages obtained from Table 2, allow us to compute an estimation of the viscoelastic effect, independent of the starting elastic model for the PSBVs, which is being used in the FEHM. 3. The improvement of injury metrics used to assess restraint systems in vehicles or the design of other preventive elements against head trauma. Currently, the estimation is often done by the injury metric called "relative motion damage measure" (RMDM) [40,41], used to predict the probability of a SDH due to the failure of PSBVs [42,43]. However, that metric was developed based on obsolete data [26], and the data from this study can be used to update that injury metric. Finally, some current limitations of the present study that could be improved in future work: Bioengineering 2021, 8, 145 12 of 15 1. The sample used is consistent but small. So, the effects of age, gender, or other anthropometric characteristics on the viscoelastic mechanical properties could not be determined. 2. In addition, a QLVE model has been used in which the relaxation function is separable, in the sense of [13]. Given the low strain levels used for the tests (since care was taken not to cause irreversible damage to the specimen from one test stage to the next), no effects of non-separability were found. However, further work could build a somewhat more general model on that basis. In any case, the proposed model is a first approximation that even explains the data that were not used for the fits (see Section 3.3). 3. Moreover, a further extension of this work would be to examine whether the re- laxation curves could be modeled, using the Prony series of stretched exponential relaxation [44]. This could lead to series with fewer and/or more accurate terms, although it is not clear if this is the case. Further work is needed in order to determine whether the use of stretched exponentials or a more general non-separable viscoelastic model would provide better models. 5. Conclusions The results reported in the current study are the first complete characterization of the viscoelastic response of PSBVs. Relaxation times and coefficients of the Prony series are computed allowing the use of these values in combination with any hyper-elastic model for PSBVs. The model explains accurately all the load cases considered. In addition, the fitted constitutive parameters can be used for computational FEHMs; these models and the obtained values could be used to reconstruct the circumstances of head trauma and for assessing the causes of TBI. Finally, the results of this study are highly relevant to assess acute SDH due to the mechanical failure of PSBVs. Author Contributions: Conceptualization, D.S.-M., S.G.-V. and C.A.-D.; methodology, formal analy- sis, project administration, supervision, validation, and visualization, D.S.-M. and S.G.-V.; software, data curation, S.G.-V. and J.L.; investigation, D.S.-M., S.G.-V., J.L. and I.G.; funding acquisition, C.A.-D. and D.S.-M.; resources I.G., J.V.-A., J.L. and M.C.R.-S.; writing—original draft preparation, D.S.-M. and S.G.-V.; writing—review and editing, I.G. and J.L. All authors have read and agreed to the published version of the manuscript. Funding: This research was funded by Fundació La Marató de TV3 (Barcelona, No. Grant 289/C/2017, 201704.30). The research was developed within the studies of the research project “Age Effects on Traumatic Brain Injuries”. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: Due to the restrictions imposed by the legal collaboration agreement between IMLFC and UPC, the data are not publicly available, although they can be accessed upon non- anonymous request to the authors of the study and under the conditions set by the aforementioned legal collaboration agreement. Conflicts of Interest: The authors declare no conflict of interest. Abbreviations The following abbreviations are used in this manuscript: PSBV(s) Parasagittal Bridging Vein(s) FEHM Finite Element Head Model QLVE Quasi-Linear Viscoelastic SDH Subdural Hematoma TBI Traumatic Brain Injury VC Viscoelastic Contribution Bioengineering 2021, 8, 145 13 of 15 Appendix A. Prediction of the USF In this appendix. we deduce what behavior is to be expected from USF(n, T ), defined in the Equation (9), and we justify mathematically why we expect a behavior given by Formula (10). The displacement δt in the cyclic test is periodic, therefore, it can be expressed by a trigonometric Fourier series: ∞ ∞ " ! !# a0 2mπ 2mπ δt = ∆δ + ∑ am cos t + ∑ bm sin t (A1) 2 m =1 T m =1 T where T = 0.55 s is the period of the cyclic stretch. The stretch λt , its corresponding rate λ̇t , the strain ε x and the strain rate ε̇ x can be calculated with the following formulas: δt δ̇t λ2t − 1 λt = 1 + , λ̇t = , εx = , ε̇ x = λt λ̇t (A2) `0 `0 2 To estimate the viscoelastic effect we consider a finite sequence of integrals of the form Z tn   (k) − tnτ−τ δ̇τ δτ Jn := e k 1+ dτ (A3) 0 `0 `0 where tn = (n − 1 + γ) T are the times for which successive force peaks occur: ∞ −(n−1+γ) τT am P1,m (θk ) + bm P2,m (θk ) + [ am P3,m (θk ) + bm P4,m ]e k ∑ (k) Jn ( T ) = (A4) m =1 (1 + 4m2 θk2 )(1 + 16m2 θk2 ) where θk = πT/τk has been introduced for brevity, and further we have Pj,m (θk ) as easily computable fourth-degree polynomials, whose coefficients depend on cos(mγπ ) and sin(mγπ ). The above formula allows us to compute successive force peaks and, thus, N ∑ gk J1 (k) 1+ (T ) k =1 USF(n, T ) = (A5) N ∑ (k) 1+ gk Jn ( T ) k =1 where the gk are the same coefficients that appear in the Equation (6) and j is the index for which | T − 1, 5τj | is minimum. It can be seen that the value of USF approaches a maximum when the θ j ∈ (4, 5) (the exact position of the maximum depends on the gk , although very little). Taking the Equation (A4) into consideration, the expression USF(n, T ) would have the form 1 + ∑kN=1 gk [c(θk ) + d(θk )] 1 + ∑kN=1 g̃k USF(n, T ) = = (A6) −n τT −n τT 1 + ∑kN=1 gk [c(θk ) + d(θk )e k ] 1 + ∑kN=1 g̃k e k where c(θk ) and d(θk ) are functions which can be computed from Equation (A4), and where the g̃k are given by gk d ( θ k ) g̃k = 1 + ∑kN=1 gk c(θk ) The expression (A6) is exact, but inconvenient, in fact, numerically it can be seen that the summand for which it is satisfied that | T − 1.5τj | is minimal, so a reasonable approximation can be obtained using only one characteristic time: ( j) 1 + g j J1 ( T ) 1 + g̃ j USF(n, T ) ≈ ( j) =   (A7) 1 + g j Jn ( T ) 1 + g̃ j exp −n τT j Bioengineering 2021, 8, 145 14 of 15 References 1. Dewan, M.C.; Rattani, A.; Gupta, S.; Baticulon, R.E.; Hung, Y.C.; Punchak, M.; Agrawal, A.; Adeleye, A.O.; Shrime, M.G.; Rubiano, A.M.; et al. Estimating the global incidence of traumatic brain injury. J. Neurosurg. 2018, 130, 1080–1097. [CrossRef] 2. Lee, M.C.; Haut, R.C. Insensitivity of tensile failure properties of human bridging veins to strain rate: Implications in biomechanics of subdural hematoma. J. Biomech. 1989, 22, 537–542. [CrossRef] 3. Sánchez-Molina, D.; García-Vilana, S.; Llumà, J.; Galtés, I.; Velázquez-Ameijide, J.; Rebollo-Soria, M.C.; Arregui-Dalmases, C. Mechanical behaviour of blood vessels during a Traumatic Brain Injury: Elastic and viscoelastic contributions. Biology 2021, 10, 831. [CrossRef] 4. Delye, H.; Goffin, J.; Verschueren, P.; V.; er Sloten, J.; Van der Perre, G.; Alaerts, H.; Verpoest, I.; Berckmans, D. Biomechanical properties of the superior sagittal sinus-bridging vein complex (No. 2006-22-0024). Sae Tech. Pap. 2006, 50, 1–12 5. Monson, K.L.; Converse, M.I.; Manley, G.T. Cerebral blood vessel damage in traumatic brain injury. Clin. Biomech. 2019, 64, 98–113. [CrossRef] 6. Monson, K.L.; Goldsmith, W.; Barbaro, N.M.; Manley, G.T. Axial mechanical properties of fresh human cerebral blood vessels. J. Biomech. Eng. 2003, 125, 288–294. [CrossRef] [PubMed] 7. Monea, A.G.; Baeck, K.; Verbeken, E.; Verpoest, I.; V.; er Sloten, J.; Goffin, J.; Depreitere, B. The biomechanical behaviour of the bridging vein–superior sagittal sinus complex with implications for the mechanopathology of acute subdural haematoma. J. Mech. Behav. Biomed. Mater. 2014, 32, 155–165. [CrossRef] [PubMed] 8. Lee, M.C.; Haut, R.C. Strain rate effects on tensile failure properties of the common carotid artery and jugular veins of ferrets. J. Biomech. 1992, 25, 925–927. [CrossRef] 9. Famaey, N.; Cui, Z.Y.; Musigazi, G.U.; Ivens, J.; Depreitere, B.; Verbeken, E.; Vander Sloten, J. Structural and mechanical characterisation of bridging veins: A review. J. Mech. Behav. Biomed. Mater. 2015, 41, 222–240. [CrossRef] [PubMed] 10. Pasquesi, S.A.; Margulies, S.S. Failure and fatigue properties of immature human and porcine parasagittal bridging veins. Ann. Biomed. Eng. 2017, 45, 1877–1889. [CrossRef] 11. Wang, Z.; Golob, M.J.; Chesler, N.C. Viscoelastic properties of cardiovascular tissues. Viscoelastic Viscoplastic Mater. 2016, 2, 64. 12. Funk, J.R.; Hall, G.W.; Crandall, J.R.; Pilkey, W.D. Linear and quasi-linear viscoelastic characterization of ankle ligaments. J. Biomech. Eng. 2000, 122, 15–22. [CrossRef] [PubMed] 13. Davis, F.M.; De Vita, R. A nonlinear constitutive model for stress relaxation in ligaments and tendons. Ann. Biomed. Eng. 2012, 40, 2541–2550. [CrossRef] [PubMed] 14. Pascoletti, G.; Di Nardo, M.; Fragomeni, G.; Barbato, V.; Capriglione, T.; Gualtieri, R.; Talevi, R.; Catapano, G.; Zanetti, E.M. Dynamic characterization of the biomechanical behaviour of bovine ovarian cortical tissue and its short-term effect on ovarian tissue and follicles. Materials 2020, 13, 3759. [CrossRef] 15. Miller, J.D.; Nader, R. Acute subdural hematoma from bridging vein rupture: A potential mechanism for growth. J. Neurosurg. 2014, 120, 1378–1384. [CrossRef] 16. Costa, J.M.; Fernandes, F.A.; de Sousa, R.J.A. Prediction of subdural haematoma based on a detailed numerical model of the cerebral bridging veins. J. Mech. Behav. Biomed. Mater. 2020, 111, 103976. [CrossRef] 17. Fung, Y.C. Biomechanics Mechanical Properties of Living Tissues; Springer Science & Business Media: New York, NY, USA, 2013. 18. Jin, Z.H. Some notes on the linear viscoelasticity of functionally graded materials. Math. Mech. Solids 2006, 11, 216–224. [CrossRef] 19. Soussou, J.E.; Moavenzadeh, F.; Gradowczyk, M.H. Application of Prony series to linear viscoelasticity. Trans. Soc. Rheol. 1970, 14, 573–584. [CrossRef] 20. Tzikang, C. Determining a Prony Series for a Viscoelastic Material from Time Varying Strain Data; NASA Center for Aero Space Information (CASI): Hampton, VA, USA, 2000. 21. Tapia-Romero, M.A.; Dehonor-Gómez, M.; Lugo-Uribe, L.E. Prony series calculation for viscoelastic behavior modeling of structural adhesives from DMA data. Ing. Investig. Tecnol. 2020, 21, 1–10. [CrossRef] 22. Sun, Z.; Gepner, B.D.; Lee, S.H.; Rigby, J.; Cottler, P.S.; Hallman, J.J.; Kerrigan, J.R. Multidirectional mechanical properties and constitutive modeling of human adipose tissue under dynamic loading. Acta Biomater. 2021, 129, 188–198. [CrossRef] 23. Braess, D. Finite Elements: Theory, Fast Solvers, and Applications in Solid Mechanics; Cambridge University Press: Cambridge, UK, 2007. 24. Iriarte, Y.A.; Varela, H.; Gómez, H.J.; Gómez, H.W. A Gamma-Type distribution with applications. Symmetry 2020, 12, 870. [CrossRef] 25. Al-A’asam, J.A. Deriving the composite Simpson rule by using Bernstein polynomials for solving Volterra integral equations. Baghdad Sci. J. 2014, 11, 1274–1282. 26. Löwenhielm, P. Dynamic properties of the parasagittal bridging beins. Z. FüR Rechtsmed. 1974, 74, 55–62. [CrossRef] 27. Sánchez-Molina, D.; García-Vilana, S.; Velázquez-Ameijide, J.; Arregui-Dalmases, C. Probabilistic assessment for clavicle fracture under compression loading: Rate-dependent behavior. Biomed. Eng. Appl. Basis Commun. 2020, 32, 2050040. [CrossRef] 28. Gefen, A.; Margulies, S.S. Are in vivo and in situ brain tissues mechanically similar? J. Biomech. 2004, 37, 1339–1352. [CrossRef] 29. Komatsu, K.; Sanctuary, C.; Shibata, T.; Shimada, A.; Botsis, J. Stress-relaxation and microscopic dynamics of rabbit periodontal ligament. J. Biomech. 2007, 40, 634–644. [CrossRef] 30. Zhang, K.; Siegmund, T.; Chan, R.W. Modeling of the transient responses of the vocal fold lamina propria. J. Mech. Behav. Biomed. Mater. 2009, 2, 93–104. [CrossRef] [PubMed] Bioengineering 2021, 8, 145 15 of 15 31. Coats, B.; Eucker, S.A.; Sullivan, S.; Margulies, S.S. Finite element model predictions of intracranial hemorrhage from non-impact, rapid head rotations in the piglet. Int. J. Dev. Neurosci. 2012, 30, 191–200. [CrossRef] [PubMed] 32. Raul, J.S.; Roth, S.; Ludes, B.; Willinger, R. Influence of the benign enlargement of the subarachnoid space on the bridging veins strain during a shaking event: A finite element study. Int. J. Leg. Med. 2008, 122, 337–340. [CrossRef] 33. Kleiven, S. Predictors for traumatic brain injuries evaluated through accident reconstructions (No. 2007-22-0003). SAE Tech. Paper 2007, 51, 1–35. 34. Doorly, M.C. Investigations into Head Injury Criteria Using Numerical Reconstruction of Real Life Accident Cases. Ph.D. Thesis, University College Dublin, Dublin, Ireland, 2007. 35. Yan, W.; Pangestu, O.D. A modified human head model for the study of impact head injury. Comput. Methods Biomech. Biomed. Eng. 2011, 14, 1049–1057. [CrossRef] [PubMed] 36. Viano, D.C.; Casson, I.R.; Pellman, E.J.; Zhang, L.; King, A.I.; Yang, K.H. Concussion in professional football: Brain responses by finite element analysis: Part 9. Neurosurgery 2005, 57, 891–916. [CrossRef] 37. Zoghi-Moghadam, M.; Sadegh, A.M. Global/local head models to analyse cerebral blood vessel rupture leading to ASDH and SAH. Comput. Methods Biomech. Biomed. Eng. 2009, 12, 1–12. [CrossRef] 38. Fernandes, F.A.; Tchepel, D.; de Sousa, R.J.A.; Ptak, M. Development and validation of a new finite element human head model: Yet another head model (YEAHM). Eng. Comput. 2018, 35, 477–496. [CrossRef] 39. Mao, H.; Zhang, L.; Jiang, B.; Genthikatti, V.V.; Jin, X.; Zhu, F.; Yang, K.H. Development of a finite element human head model partially validated with thirty five experimental cases. J. Biomech. Eng. 2013, 135, 111002. [CrossRef] [PubMed] 40. Sanchez-Molina, D.; Velazquez-Ameijide, J.; Arregui-Dalmases, C.; Crandall, J.R.; Untaroiu, C.D. Minimization of analytical injury metrics for head impact injuries. Traffic Inj. Prev. 2012, 13, 278–285. [CrossRef] 41. Sánchez-Molina, D.; Arregui-Dalmases, C.; Velázquez-Ameijide, J.; Angelini, M.; Kerrigan, J.; Crandall, J. Traumatic brain injury in pedestrian–vehicle collisions: Convexity and suitability of some functionals used as injury metrics. Comput. Methods Programs Biomed. 2016, 136, 55–64. [CrossRef] [PubMed] 42. Takhounts, E.G.; Ridella, S.A.; Hasija, V.; Tannous, R.E.; Campbell, J.Q.; Malone, D.; Duma, S. Investigation of traumatic brain injuries using the next generation of simulated injury monitor (SIMon) finite element head model (No. 2008-22-0001). SAE Tech. Paper 2008, 52, 1. 43. Fernandes, F.A.; Sousa, R.J.A.D. Head injury predictors in sports trauma–a state-of-the-art review. Proc. Inst. Mech. Eng. Part J. Eng. Med. 2015, 229, 592–608. [CrossRef] 44. Mauro, J.C.; Mauro, Y.Z. On the Prony series representation of stretched exponential relaxation. Phys. Stat. Mech. Its Appl. 2018, 506, 75–87. [CrossRef]