A Data-Driven Memory-Dependent Modeling Framework for Anomalous Rheology: Application to Urinary Bladder Tissue

We introduce a data-driven fractional modeling framework aimed at complex materials, and particularly bio-tissues. From multi-step relaxation experiments of distinct anatomical locations of porcine urinary bladder, we identify an anomalous relaxation character, with two power-law-like behaviors for short/long long times, and nonlinearity for strains greater than 25%. The first component of our framework is an existence study, to determine admissible fractional viscoelastic models that qualitatively describe the linear relaxation behavior. After the linear viscoelastic model is selected, the second stage adds large-strain effects to the framework through a fractional quasi-linear viscoelastic approach, given by a multiplicative kernel decomposition of the selected relaxation function and a nonlinear elastic response for the bio-tissue of interest. From single-relaxation data of the urinary bladder, a fractional Maxwell model captures both short/long-term behaviors with two fractional orders, being the most suitable model for small strains at the first stage. For the second stage, multi-step relaxation data under large strains were employed to calibrate a four-parameter fractional quasi-linear viscoelastic model, that combines a Scott-Blair relaxation function and an exponential instantaneous stress response, to describe the elastin/collagen phases of bladder rheology. Our obtained results demonstrate that the employed fractional quasi-linear model, with a single fractional order in the range $\alpha = 0.25 - 0.30$, is suitable for the porcine urinary bladder, producing errors below 2% without need for recalibration over subsequent applied strains. We conclude that fractional models are attractive tools to capture the bladder tissue behavior under small-to-large strains and multiple time scales, therefore being potential alternatives to describe multiple stages of bladder functionality.

1. Introduction. Bio-tissues are complex and multi-functional materials, optimized for their specific host organisms, and constrained by limited set of building blocks and available resources [26]. While the mechanical behavior of a number of standard engineering materials is quite well-understood, there is still a significant effort towards bio-materials, where microstructure heterogeneities, randomness and small scale physical mechanisms lead to non-standard and at times counter-intuitive responses. Power-law viscoelastic rheology is a complex response observed in many bio-tissues such as arteries [15], cartilage [38], lungs [59], smooth muscle [19], liver and kidneys [49], among other classes of materials. These power-law materials, also termed anomalous, exhibit one or more power-law scalings for creep/relaxation in the form J(t) ∝ t β and G(t) ∝ t −β across multiple time-scales. Similar anomalous behav-iors are also present for dynamic storage/dissipation in the frequency domain [40,64]. The origin of this power-law behavior at the continuum level is linked to (non-Fickian) sub-diffusive processes [41] in the corresponding fractal-like micro-structures [5].
The aforementioned anomalous non-exponential behavior usually requires a significant number of material parameters when employing standard viscoelastic models. These consist of mechanical arrangements of linear springs and Newtonian dashpots, which induces a finite number of relaxation modes, which may lack predictability when performing outside the experimental time scales [27,38]. In this regard, fractional models become attractive alternatives, since their integro-differential operators naturally utilize power-law convolution kernels, coding self-similar microstructural features in a reduced-order mathematical language with smaller parameter spaces. Therefore, they have been employed as compact and predictive models for a number of anomalous systems, such as biological materials [38,15,59,19,49], fluid turbulence [53,54,1], and instabilities [2]. We particularly note that such predictability has been shown to extend across different experiments (relaxation/creep) in certain cases [27]. Additionally, calibrating experimental data with a set of existing rheological models leads to a material model selection problem, which is inherently ill-conditioned, since multiple models can pragmatically yield similar errors when confronted to experiment. In this work, we attempt to reduce this implicit ill-posedness by introducing fractional-order models as attractive alternatives to their integer-order counterparts, and employed to urinary bladder (UB) tissue modeling. Our fractional modeling framework aims to obtain compact mathematical models with a reduced number of material model parameters, while introducing a minimal, but sufficient number of fractional rheological elements that capture the qualitative response of multiple power laws and minimizes the errors, also rigorously taking into account the corresponding power-law memory effects.
The lower urinary tract, and especially the UB, is a highly dynamic organ system. To ensure its proper function, the bladder needs to be able to significantly increase in size while maintaining a low internal pressure, and this ability is dictated by the mechanics of the bladder wall. Specifically, during filling, the bladder tissue must leverage its viscoelastic characteristics to accommodate for large deformations without resulting in significant increase of luminal pressure. When this behavior is compromised due to disease, the resulting increase in pressure might generate a highpressure urine reflux from the bladder to the kidneys, resulting in renal failure [6,37]. To increase the complexity of the organ mechanics, the characteristics of the bladder differ between different anatomical locations (i.e., dorsal, ventral, lateral, lower-body, trigone) [34,42] and orientations (i.e., longitudinal/apex-to-base and circumferential/transverse) [34,42,11,22,12,68,48,76,47]. To describe the mechanical behavior of bladder tissue, both hyperelastic [12,76,47,22,33,16,70,67,51,17,75,24] and viscoelastic [68,48,13,66,14,69,23,71,4,60,32,3,72,73,43,44] models have been used in the literature. However, due to the differences in mechanical testing protocols as well as modeling, most of the results cannot be compared with one another and often results in contradicting conclusions. While several pathologies of the lower urinary tract are associated with dramatic changes of the mechanical behavior of the bladder wall [65], still much is unknown about the mechanisms that affect this organ, not just in diseased states but in healthy as well. In this study, we focus on the healthy behavior of the porcine urinary bladder, which a present work suggested is a good model for the human urinary bladder.
Although fractional linear viscoelasticity has been succesfully employed in a number of biomechanical applications, additional modeling considerations are necessary when dealing with other material nonlinearities, such as large strains. This would imply that the material relaxation behavior depends on both time and applied strains, which requires additional modeling considerations. In the UB case, Korossis et al. [34] reports the small-strain regime to be drive by elastin, while a stiffer response for large strains was driven by collagen. Jokandan et al. [28] characterized the quasistatic stress-strain response of the UB as exponential-like. To address such behavior, a practical and well-known class of models is the quasi-linear-viscoelastic (QLV) theory by Fung [21], which considers a multiplicative coupling between a linear viscoelastic relaxation and a nonlinear elasticity term. While the original formulation utilizes a finite spectrum of relaxation times, fractional extensions of the QLV theory have been developed and employed for the response of bio-tissues [20,15]. Doehring et al. [20] employed a Mittag-Leffler-type reduced relaxation function that captured the short and long term behaviors of aortic valve cusp. Craiem et al. [15] developed a fractional Kelvin-Voigt-type reduced relaxation function with an exponential instantaneous response, which was successfully applied to the nonlinear viscoelasticity of arterial walls. Regarding the nonlinear viscoelastic modeling of the UB, Natali et al. [48], developed an anisotropic, visco-hyperelastic model, which was validated through a uniaxial experiment under cyclic loads. Nagatomi et al. [43] studied the nonlinear behavior of rat bladders, by calibrating a two-dimensional QLV model to bi-axial relaxation data from bladders of subjects that were healthy and with spinal cord injury. Their findings reported a need for new models to account for both normal and pathological states, due to tissue remodeling.
To the authors' best understanding, although existing studies have addressed the nonlinear viscoelastic of the UB for different subjects and loading conditions, there are no studies in the literature leveraging the use of fractional viscoelastic models to model potential, emerging power-law behaviors. In this work we develop a datadriven fractional modeling framework for linear and quasi-linear viscoelasticity to account for both anomalous power-law relaxation and large strains of bio-tissues. We validate the developed framework for the first time in the uniaxial relaxation of porcine urinary bladder tissue for a wide range of applied strains. The characteristics of our experimental procedure are: • We obtain the porcine UB uniaxial relaxation data from small-to-large strains of five distinct anatomical locations. • Our relaxation experiments are performed under increasingly larger strains, without intermediate unloading steps or tissue preconditioning. • The mechanical response of the UB indicates nonlinear viscoelastic behavior with power-law relaxation, characterizing an anomalous, non-exponential behavior. Given the anomalous nonlinear response of the UB tissue, we develop our two-stage anomalous modeling framework as follows: • In the first stage, we develop an existence study that considers a set of linear fractional building block models (Scott-Blair, fractional Kelvin-Voigt, fractional Maxwell), which are selected according to the multi-power-law nature of the relaxation data and calibration errors at the linear viscoelastic regime. • In the second stage, we account for the large strain behavior of the corresponding tissue, by employing a fractional quasi-linear viscoelastic (FQLV). The goal is to extend the quality of power-law relaxation (due to the material's fractal microstructure) to the large strain regime of the tissue of interest. We employed the aforementioned two-stage formulation to the UB experimental data, and obtained the following main findings; • All candidate fractional linear viscoelastic models provide sufficiently accurate fits for single-relaxation steps under smaller strain levels, where the two fractional order Maxwell model is the most suited for the UB data. • The employed four-parameter FQLV model with a reduced Scott-Blair relaxation function and exponential instantaneous stress response was successful over five consecutive relaxation steps, with root mean squared errors below 2%, and without the need of model recalibration between applied step strains. • The lower-range of obtained fractional orders is around α = 0.17−0.30, which is compatible with the observed long-time slopes of the UB relaxation data.
Small α values have been suggested to indicate strong fractality in bio-tissue microstructures such as collagen fibers [20]. The rest of the paper is organized as follows. In Section 2 we present the problem setup and methodology, with the uniaxial UB stress relaxation experiments and our proposed fractional modeling framework for biotissues, comprised of linear and quasilinear fractional models. In Section 3, we present our obtained linear viscoelasticity results for the UB relaxation under the first strain step, and the fractional quasi-linear viscoelastic model for all consecutive strain steps, followed by the discussions and potential improvements of the work. Finally, in Section 4, we provide the concluding remarks of this work and future directions.  30,45,45,45, and 45 minutes, respectively. Besides the strain inputs, the force denoted by F data (t) is measured by a 10 [lb] load cell throughout the duration of the test. The cross-sectional area A data (t) of each sample is calculated by taking top and side view pictures that are converted to binary images which are processed in MATLAB ® , respectively estimating the base b(t) and height h(t) dimensions of the samples after the application of each consecutive step strain ε i . The updated crosssectional area is assumed to remain constant throughout the relaxation at each strain level, and is evaluated as A data = b(t)h(t). Given force and area time-series, the true strain is evaluated as σ data (t) = F data (t)/A data (t).
Once the stress σ data (t) and strain ε data (t) time-series are obtained, we filter the data through a moving average filter with a time-window of thirty neighbor data points. Figure 2.2 illustrates the relaxation curves for all samples in linear and log-log scales. We observe a characteristic power-law scaling for long-time behavior, which is evident in Fig.2.2 (b). As will be shown later through fractional model fits, the relatively low scaling coefficient β indicates an anomalous behavior of predominantly elastic nature, with a plateau with low decay rates σ ∼ t −β at larger time-scales (i.e., t > 400 [s]). We also note that the trigone and lower body specimens yielded higher stress levels, particularly at very high strains, while the dorsal specimen yielded lower overall values. This is in accordance to stress-strain results obtained by Korossis et al. [34], that indicated statistically significant, higher collagen phase slopes for the lower body and trigone regions. We analyze the presence of strain dependency on the UB relaxation behavior, in order to determine if the viscoelasticity is of linear or  nonlinear nature. Therefore, we employ the definition of linear relaxation modulus G(t) [P a], applied for each fixed strain application from experimental data [25]: with i = 0, 1, . . . , N steps . In addition to the strain dependency on the relaxation behavior, the above definition also allows us to identify the presence of additional power-laws or stretched exponential behaviors in the tissue response. We observe that although the relaxation moduli data for each sample is almost linear for ε i = 0.25, 0.50, since the curves approximately overlap (except for the trigone sample), the behavior of G data is clearly time-and strain-dependent. Furthermore, the degree of nonlinearity is more pronounced for the lower body and trigone samples, and less pronounced for the dorsal sample. Interestingly, we notice two limiting powerlaw behaviors for short and long times. The short time behavior appears to have a stretched exponential nature, with a limiting power law of smaller magnitude β 1 , while the long time behavior is associated with a power-law of larger magnitude β 2 , see Fig.2.3(f). Nevertheless, the relaxation behavior clearly transitions from a slower regime to a faster regime, even though the latter still presents a far-from-equilibrum response (no equilibrium glass state) within the experimental time scales. Larger standard deviations for the long time power-law were observed for the trigone region due to its distinct response for ε 0 = 0.25. We remark that this analysis is just performed to infer the nonlinear relaxation quality the data, and we do not intend to construct a master curve for each UB sample.

A Fractional Viscoelastic
Modeling Framework for Anomalous Tissue Rheology. We develop a two-stage framework to assess the qualitative linear/nonlinear relaxation behavior of anomalous materials and identify the most feasible fractional viscoelastic models, in order to alleviate the inherent ill-posedness of model selection problems. In the first stage, we develop an existence study to identify the admissible set of anomalous constitutive laws that satisfy the quality of the experimental linear relaxation behavior, while shedding light on the corresponding microstructural constituents associated to anomalous behavior. In the second stage, we incorporate large strain effects to the relaxation quality, based on the nonlinear nature of the tissue of interest.

First stage: An existence study of fractional linear viscoelastic models.
Starting with the Scott-Blair (SB) model as the fundamental building block, we construct building block models through parallel and serial combinations to obtain the fractional Kelvin-Voigt (FKV) and fractional Maxwell (FM) models. In our approach, we take into account the anomalous qualities present in the experimental relaxation data and compare them with each of the candidate building block models. Each of the models exhibit distinguished material complexities, such as distinct asymptotic behaviors of relaxation G(t), multiple power-law regimes, slower/faster relaxation at the asymptotic stages [55,39] and presence of material nonlinearities. In the last part of the existence study we classify the candidate models according to their anomalous response, which together with obtained fitting errors and number of material parameters, constitutie the criteria for the model selection procedure. Given the experimental data presented in Section 2.1, we focus on the first relaxation step (ε 0 = 0.25) of Figs.2.2 and 2.3 for the first stage of our data-driven framework. Our objective is to demonstrate how fractional viscoelastic models are able to capture the UB relaxation with simplistic mechanical arrangements and a small number of material parameters.
The rheological building block for our framework is the fractional SB viscoelastic element, which compactly represents an anomalous viscoelastic constitutive law connecting the stresses and strains: with t > 0, ε(0) = 0, and constant fractional order in the range 0 < α < 1, which provides a material interpolation between Hookean (α → 0) and Newtonian (α → 1) elements. The operator C 0 D α t (·) represents the time-fractional Caputo derivative given by: where Γ(·) represents the Euler-gamma function [39] andu = du/dt. The pair (α, E) uniquely represents the SB constants, where the pseudo-constant E [P a.s α ] compactly describes textural properties, such as the firmness of the material [8,27]. In this sense E is interpreted as describing a snapshot of a non-equilibrium dynamic process instead of an equilibrium state. The corresponding rheological symbol for the SB model represents a fractal-like arrangement of springs and dashpots [56,40], which we interpret as a compact, upscaled representation of a fractal-like microstructure. Regarding the thermodynamic admissibility of the SB element and more complex models (i.e., plasticity and damage) involving it, we refer the reader to Lion [36] for the SB model, and Suzuki et al. [63]. The relaxation function G SB (t) [P a] for the SB model is given by the following single, inverse power-law form: which is the convolution kernel of the integro-differential form in (2.3), with a modulating pseudo-constant E for fixed α. which is scale-free, i.e., a single power-law is present for all t > 0. We note that although this relaxation response may seem to be oversimplified, it provides a flexible constitutive interpolation able to, at the very least, take into account the long-term anomalous dynamics of materials, such as the power-law β 2 in Fig.2.3. This also allows the SB element to capture, in certain time-scales, power-law behaviors induced by predominantly elastic microstructures, such as collagen networks [20] with small α-values. We also note that single collagen fibril relaxation data [57] has been calibrated to single-order fractional Kelvin-Zener models (where a single SB element is combined with Hookean springs), with reported values in the range α = 0.12−0.21 [9]. Similarly, the macroscopic response of ultrasoft tissues such as the brain, which has a more fluid-like character with significant microscopic inter-cell body displacements in the deformation process [52], has also been captured using SB elements combined with Hookean springs [31,18]. Kohandel et al. [31] calibrated a single-order dynamic modulus of a fractional Kelvin-Zener model to bovine brain tissue, obtaining a fractional order α = 0.6. Through relaxation experiments, Davis et al. [18] calibrated a similar model to bovine brain paerenchima also reporting a high value of α = 0.64.
We utilize the SB model as our rheological building block, and define a set of "building block models", which introduce a higher degree of material complexity through multiple power-law behaviors for relaxation and therefore distinct anomalous regimes for small and large time-scales. This multi-fractal type of rheology is characteristic of cells [58] and biological tissues [74], due to their complex, hierarchical and heterogeneous microstructure. Here we consider the two simplest canonical combinations of SB elements. Through a parallel combination, we obtain the fractional Kelvin-Voigt (FKV) model, which has the following stress-strain relationship [56]: , with t > 0, ε(0) = 0, fractional orders 0 < α 1 , α 2 < 1 and associated pseudo-constants E 1 [P a.s α1 ] and E 2 [P a.s α2 ]. The corresponding relaxation modulus G(t), is also an additive form of relaxation involving two SB elements: We notice a response characterized by two power-law regimes, with a transition from faster to slower relaxation slopes. The asymptotic responses for small and large time-scales are given by G FKV ∼ t −α2 as t → 0 and G FKV ∼ t −α1 as t → ∞. We note that this quality allows the FKV model to describe materials that reach an equilibrium behavior for large times when α 1 → 0, which is intuitive from the mechanistic standpoint as one of the SB elements becomes a Hookean spring. Regarding applications of the FKV model to bio-tissue constituents, Bonfanti et al. [9] has found the combination of low-high fractional orders to be proper for bovine trachea smooth muscle cells, recovering α 1 = 0.1 and α 2 = 0.78. Finally, through a serial combination of SB elements, we obtain the fractional Maxwell (FM) model [27], given by: with t > 0, fractional orders 0 < α 1 < α 2 < 1, the additional constraint 0 < α 2 −α 1 < 1, and two sets of initial conditions for strains ε(0) = 0, and stresses σ(0) = 0. We note that in the case of non-homogeneous ICs, one requires compatibility conditions [39] between stresses and strains at t = 0. The corresponding relaxation function for this building block model assumes a more complex, Miller-Ross form [27]: where E a,b (z) denotes the two-parameter Mittag-Leffler function, defined as [39]: Interestingly, the presence of a Mittag-Leffler function in (2.8) leads to a stretched exponential relaxation for smaller times and a power-law behavior for longer times, as illustrated in Fig.2.4. We also observe that the limit cases are given by G FM ∼ t −α1 as t → 0 and G FM ∼ t −α2 as t → ∞, indicating that the FM model provides a behavior transitioning from slower-to-faster relaxation. Furthermore, when α 2 → 1, the FM model presents a fluid-like behavior for long times [55], and therefore allowing a fast relaxation, similar to its integer-order counterpart. We refer the reader to the recent work by Bonfanti et al. [9] for a number of applications of the aforementioned models. We notice that both FKV and FM models are able to recover the SB element with a convenient set of pseudo-constants. Furthermore, from Fig.2.4, we observe that, under a variable order α(t) setting, the SB relaxation function (2.4) is able to recover both constant-order FKV and FM models, respectively, through parametric decreasing/increasing variable order functions. However, since variable-order models may potentially increase the material parameter space and computational cost, we restrict our discussion to constant order models. Furthermore, we also outline more complex building block models that yield more flexible responses, including three to four fractional orders, such as the fractional Kelvin-Zener (FKZ), fractional Poynting-Thomson (FPT), and fractional Burgers (FB), which in turn are able to recover the FKV and FM models. We refer the reader to the works [56,9] for more details on such models.

Second stage:
Fractional quasi-linear viscoelastic modeling. The presented models in Section 2.2.1 provide candidates for power-law relaxation functions that describe the anomalous linear viscoelasticity of biotissues, however, in such materials the stress-strain relationship may becomes nonlinear as fibers in collagen networks transition from entangled to aligned with the applied load direction. Therefore, the viscoelastic behavior itself becomes nonlinear and the relaxation function has an intrinsic dependency on the strain levels, as observed in Fig.2.3 under successive large step-strain applications. To incorporate this additional effect to our modeling framework, we follow [21,15], and employ the following quasi-linear, fractional viscoelastic model (FQLV): where the convolution kernel is given by a multiplicative decomposition of a reduced relaxation function G(t) which may be selected from the first stage of the framework, and an instantaneous, nonlinear elastic tangent response with stress σ e . In the work by Craiem et al. [15], the reduced relaxation function has a fractional Kelvin-Voigtlike form with one of the SB elements replaced with a Hookean element. Here, we assume a simpler rheology and adopt a Scott-Blair-like reduced relaxation in the form: with the pseudo-constant E with units [s α ], since the elastic strains σ e (ε) have units of [P a]. We adopt the same, two-parameter, exponential nonlinear elastic part as in [15]: with A having units of [P a] and B being a non-dimensional tuning parameter for the degree of nonlinearity induced by applied strains. Substituting Equations 2.11 and 2.12 into Eq.2.10, we obtain: which differs slightly from the linear SB model (2.2) in the sense that an additional exponential factor multiplies the function being convoluted. We remark that reduced relaxation forms for the FKV and FM can also be employed in the FQLV framework. For the FKV model, it would lead to a multi-term convolution of the same nature as (2.13) with fractional orders α 1 and α 2 . Therefore, the same numerical methods employed in Section 2.2.3 could be employed. However, when dealing with a FM reduced relaxation function, a specialized type of quadrature for a Miller-Ross-type of kernel as (2.8) would be required, potentially with impractical computations of Mittag-Leffler functions for a large number of time-steps.

Numerical discretizations.
We discretize the fractional Caputo derivatives in Equations 2.2-2.7 through an implicit L1 finite-difference scheme [35]. Therefore, we consider a uniform time-grid with N time-steps of size ∆t, such that t n = n∆t, with n = 0, 1, . . . , N . We remark that although the equations for each of the building block models could be discretized utilizing fast schemes and approaches that treat initial time singularities (see [78,79] and references therein), the number of utilized data-points is not large, with N data ≈ 3 000 for the first relaxation and N data ≈ 25 000 for all steps. Also, we avoid the singularity nearby t ≈ 0 since the first relaxation step is applied at approximately 80 seconds. Nevertheless, the non-smooth nature of the loading would degenerate most of the existing numerical methods for FDEs, and we found that the employed method in this work with the ∆t described in Section 3 is sufficient for the accuracy to reach the plateau of the experimental data, such that model error is dominant.
In the following, we present the discretized forms for each of our employed linear fractional viscoelastic models. We start with the stresses for the SB model evaluated at t = t n+1 : (2.14) σ For the FKV model, we obtain: with discretization constants Finally, for the FM model, we have: with discretization constants The history terms H ν u with fractional order ν in the above equations are given by the following form: with weights b ν j := (j + 1) 1−ν − j 1−ν . The discretization for the FQLV model (2.13) employed in this work is shown in [61], which is a straightforward, fully-implicit, L1 finite-difference approach with a trapezoidal rule employed on the additional exponential factor. Therefore, the discretized stresses for the FQLV model are given by: with constant C 1 = EAB/(∆t α Γ(2 − α)). The discretized history load in this case is given by: with discretization weights b α k = (k + 1) 1−α − k 1−α and ε i+ 1 2 = (ε i + ε i+1 )/2. The presented discretization has an accuracy of O(∆t 2−α ), and we refer the reader to [61] for simulations of numerical convergence.

Model optimization.
We perform the model fits through a particle-swarm optimization (PSO) algorithm [29], which was implemented in MATLAB. The adopted PSO parameters are a population N pop = 30 and N it = 1000 iterations for the linear cases and N it = 100 iterations for the nonlinear cases. For the linear viscoelastic fits, we set the initial material pseudo-parameter in the 0 ≤ E ≤ 10 8 [P a.s α ], and the fractional orders are constrained in the 0.0001 ≤ α ≤ 0.9999 range, to ensure that the employed fractional models are able to recover simpler fractional counterparts and also standard rheological elements, if required by the experimental data. For nonlinear cases, we estimate ranges for parameters A and B of the FQLV model by fitting the instantaneous stress response (2.12) to each stress peak in every step-strain application of Fig.2.2. From this preliminary estimate, we have obtained parameters in the ranges 10 4 ≤ A ≤ 10 5 [P a] and 0 ≤ B ≤ 2, which are taken as input parameter ranges for the PSO algorithm. For the relaxation parameters of the FQLV, as in [15], we note that the nature of the power-law relaxation kernel, it is nontrivial to obtain a normalized G(0 + ) = 1. Nevertheless, for the pseudo-constant we set the range 0 ≤ E ≤ 1 and for the fractional order α we employ the same range as the linear case.
Since the stresses σ data (i) and strains ε data (i) from the relaxation dataset are non-uniform in time, we perform a linear (first-order accurate) interpolation of the strains ε data (i) to an uniform grid. We then utilize the input strains and compute the global best solution for stress for every PSO iteration through (2.14), (2.15), (2.16), or (2.17). Then, we linearly interpolate the stress back to the nonuniform grid to obtain σ model . The time-step size for the uniform grid solution is set to ∆t = 0.495 [s], which is the minimum time interval between two consecutive datapoints. For verification purposes, we tested smaller step-sizes (∆t = 0.0495) and did not obtain improved results, and we note that all employed numerical discretizations for fractional models are fully-implicit. Therefore, our verification step indicates that model error dominates over discretization error for the employed time-step size. The cost function is defined as: The adopted error measures in this work are the normalized least-squares error (LSE) and root mean squared error (RMSE) between the experimental and mapped simulated stresses, which are respectively given by: Finally, all numerical simulations were run in a computer system with Intel Xeon Gold 6148 CPUs with 2.40GHz.

Results and Discussion.
3.1. Linear Viscoelasticity. Figure 3.1a illustrates the obtained fits for all bladder samples utilizing an SB model, for the first strain step (ε 0 = 0.25). We observe very good fits for most samples, especially at larger time scales, with an exception for the trigone (T) sample due to a sudden stress drop in the experimental data. The fitting quality decreases for all samples at the early relaxation dynamics (nearby the step-strain application), with the SB model underestimating the maximum values of stress peaks. The obtained fractional orders lie in the 0.2 − 0.3 range (see Table 3.1), which are similar to the observed long-term power law from the estimated experimental relaxation functions in Fig.2.3. Furthermore, the least-squared errors lie in the 2% − 11% range, while the RMS errors are within the 2% − 4% range. The higher values of the pseudo-constant and fractional order for the trigone sample are likely due to the SB model accounting for both instantaneous and long-term response over its limited set of two parameters. We also note that the FKV model pragmatically recovered the SB model in all instances, where the PSO algorithm obtained optimal values for the fractional orders that are close to the SB model. In addition, the optimal values for one pseudo-constant is either set to zero, or the summation of E 1 and E 2 recovers the value of E for the SB model. This indicates the a FKV model does not improve the bladder fit quality, and one would rather employ a SB model with half the amount of material parameters under the same error levels.     Table 3.1 indicate the presence of a predominantly elastic power-law α 1 in the 0.17 − 0.19 range, and a predominantly viscous α 2 in the 0.74 − 0.99 range. Particularly, the FM model fit for the trigone specimen indicates the recovery of a dashpot element, and thus the corresponding SB element could be replace by a Newton element. Regarding pseudo-constant values, we note that E 1 values have variations that qualitatively agree with the intensity of stress peaks, but E 2 values can vary in several orders of magnitude, which could be due to the presence of multiple local minima or the discrepancy between obtained fractional orders α 2 . In general, the dorsal and ventral samples seem to be the most anomalous, as they present both fractional order values sufficiently far from standard elements. Figure 3.2 illustrates the pointwise relative errors between the SB and FM models and the experimental data for the first strain step. We notice that the SB element has similar errors as the FM model in the 800 < t < 1200 [s] range and larger errors for longer times for most samples. For shorter time ranges, the SB model has larger errors (up to 1 order of magnitude) for all samples. This reinforces the fact that the FM model is more descriptive of both early and long-term dynamics of bladder  relaxation, as the qualitative analysis and estimated experimental relaxation moduli suggest. Furthermore, we also note that the better performance of the FM model is also attributed to better approximating the loading ramp and the peak stress preceding the relaxation behavior.
We also employed more complex fractional linear viscoelastic models, such as fractional Kelvin-Zener, Poynting-Thomson and Burgers' (see [61] for the models and their corresponding discretizations). From our employed fitting procedure, all models either recovered or had the same performance as the FM model.
3.2. Nonlinear Viscoelasticity. Figure 3.3 illustrates the obtained fits for the fractional QLV model under all consecutive strain steps, where we observe a very good agreement with the experimental data. Except for the trigone sample, all cases had higher deviations towards the final strain steps. Nevertheless, we note that the error levels are below 6% (LSE) and 2% (RMSE) for the entire dataset, under 4 material parameters, which are listed in Table 3.2. Furthermore, the obtained fractional-orders lie in the range 0.24 − 0.3 which are in accordance with the estimated power-laws in our a-priori analysis presented in Fig.2.3. Particularly, the lowest fractional order was obtained for the trigone specimen, and highest for the dorsal one. A slightly higher degree of nonlinearity is also recovered for the trigone and lower-body samples due to the larger values of B.   Figure 3.4 illustrates the pointwise relative errors for the FQLV model under all bladder samples and strain steps. We observe a similar error behavior as the SB model in Fig.3.2, particularly for larger applied strains.
3.3. Discussion. To our best understanding, this was the first work in the literature addressing fractional viscoelastic modeling to bladder tissues. From the results obtained for our building block models, the overall lower range of fractional orders obtained for all linear/nonlinear models is 0.17 − 0.3, indicating a predominantly elastic yet highly anomalous behavior with smaller decay rates at long times, i.e., the presence of far-from-equilibrium dynamics. A similar parametric range was obtained in other anomalous systems such as arterial wall relaxation [15], aortic valve tissue [20], 1D and 3D brain artery walls under fluid-structure interactions [50,77], canine and bovine liver tissue [30,10], and lung tissue [59]. As suggested by Doehring et al. [20], small α-values can be indications of strong fractality in bio-tissue microstructure such as collagen fibers, which are vastly present in the UB, and particularly with a larger network in the trigone region. The larger values of fractional orders α 2 in the 0.74 − 0.99 range obtained by the fractional Maxwell model is similar to those obtained for brain tissue relaxation [18] and human ear [46,45]. This indicates a significantly more dissipative behavior, possibly compensating the highly-anomalous behavior provided by the smaller fractional order α 1 for short time-scales, and thus better fitting the slower relaxation nearby the load application. The transitional behavior from slower-to-faster relaxation slopes observed from the UB specimens and captured by the FM model were also noticed in bio-tissues composed of weakly crosslinked collagen networks [74]. We note that although the captured fractional orders α 1 , α 2 for the FM model on the UB relaxation do not quantitatively match the slopes β 1 , β 2 for G data (t) in Fig.2.3, these fractional orders refer to the asymptotic behaviors at t → 0 and t → ∞, as illustrated in our existence study in Fig.2.4, and it is likely that relaxation experiments under a larger range of time-scales would yield a better quantitative agreement. For the purpose of our existence study, we consider a qualitative agreement and small error levels to be sufficient to select a valid candidate building block model.
The nonlinear viscoelastic behavior was well approximated by the employed FQLV model, which decomposes the relaxation kernel in a multiplicative fashion into a power-law reduced relaxation function and a tangent elastic stiffness described by an exponential elastic stress form. This allowed the nonlinear part of our existence study to capture the complex rheology of the UB with large applied strains (up to 200%) and RMS errors as low as 1%. In fact, Jokandan et al. [28] observed an exponential-like stress-strain response in quasi-static tensile testing of porcine bladder samples. Specifically, under relaxed states, the entangled configuration of collagen fibers yield a linear stress strain relationship, but a nonlinear regime with much higher stress levels is attained once the fibers align with the load direction and store most of the strain energy in the system. Korossis et al. [34] attributed the linear region to be predominantly driven by elastin, and the nonlinear phase by collagen.
Our pointwise relative error analysis reinforced the idea that the FM model is more descriptive of both slower early dynamics of the porcine bladder, but also the faster dynamics observed at longer time-scales in the estimated experimental relaxation modulus. In our error analysis for the nonlinear case, Fig.3.4 indicates a similar qualitative error behavior between the FQLV model and the linear SB element in Fig.3.2, especially towards the larger strain regime, with higher errors in the small and large times after the step-strain applications. This suggests that coupling a fractional Maxwell-type reduced relaxation function to the existing framework would very likely improve our results for the nonlinear case as it did in the linear one.
Regarding our developed framework, the existence study proved to be interesting to identify the most proper fractional linear viscoelastic model for stress relaxation, which can later inform the fractional quasi-linear viscoelastic model on the proper form of the reduced relaxation function. For the UB, we conclude that while the SB, FKV and FM models yield errors in the same order of magnitude, the FM model better captures the two power-law qualitative behavior of the data, which is fundamental for both short-and long-term predictions of tissue response. Nevertheless, the SB model provided satisfactory results for the observed experimental time-scale, and the FKV model proved to be redundant and a source of ill-posedness in a model selection framework, since it obtained the same performance as the SB model with twice as many parameters. Although we cannot guarantee that the obtained model parameters provide a global minimum for the cost function (2.19), we find our obtained fitting errors, increased number of material parameters, and diverging qualitative behaviors between the experimental data in Fig.2.3 and the FKV relaxation behavior from Fig.2.4 to be sufficient to exclude the FKV as a viable candidate for the UB. The same analysis applies to other tested models not shown here, such as fractional Kelvin-Zener, fractional Poynting-Thomson and fractional Burgers' models, which consistently recovered the FM model.
Regarding potential improvements, anisotropy of bladder samples has been observed in existing studies [42], and therefore an interesting step would be to map the variation of fractional-orders for distinct sample orientations. Based on the obtained linear and nonlinear results, another interesting aspect would be to incorporate a fractional Maxwell-type relaxation function to the fractional QLV framework, similar to Doehring et al. [20] to better describe the initial relaxation process in each strain step. However, due to the larger number of time-steps required by our dataset, an efficient numerical method would be required to handle the resulting differ-integral with a Miller-Ross relaxation kernel, or the FQLV framework would need to be developed in differential form, likely as a system of equations comprised of an FDE solving for elastic stresses and a separate equation for nonlinear elasticity. Finally, biaxial tests and models would give insight in the effects of shear stress to the tissue behavior, and confronting creep predictions with experiments would allow one to verify the consistency of the obtained parameters, as already succesfully done with other anomalous materials [27].

Conclusions.
We developed a data-driven fractional modeling framework for linear and nonlinear viscoelasticity which was validated for the first time in the uniaxial relaxation of porcine urinary bladder tissue for a wide range of applied strains. Our approach employed fractional linear and quasi-linear viscoelastic models to account for anomalous power-law relaxation and large strains. Our main findings in this study were: • Our existence study was able to relate the power-law features of specimens from five distinct bladder samples to the fractional orders in linear/quasilinear fractional models, and is an interesting step towards an automated model selection framework. • The bladder uniaxial relaxation data was obtained from consecutive and increasing step strain applications, indicating the presence of nonlinear, straindependent effects on the relaxation functions. • Our obtained data is consistent with other bladder rheology studies, indicating higher stress levels for the trigone region, and an exponential-like stressstrain relationship. • Among linear viscoelastic models employed for the first relaxation step (25% strains), the fractional Maxwell model was the most suited for all regional bladder samples, with two fractional orders, which dominate short-and longtimes. More complex linear fractional models consistently recovered the FM, and the FKV model proved to be unfeasible since it recovers a SB element. • The employed fractional quasi-linear viscoelastic model successfully captured the multi-step relaxation behavior with four material parameters and without any requirement of parameter recalibration, yielding a fractional order range α = 0.25 − 0.30 with root mean squared errors below 2%. • Fractional calculus can be a interesting modeling alternative to describe the linear/nonlinear behavior of porcine UB especially within a material model selection framework, since fractional models potentially provide a reduced number of material parameters due to the presence of multiple power-laws in relaxation. Regarding potential future steps, investigating the possibility of plastic deformations under large strains by employing quasi-linear visco-plastic effects [61,62,63] and also failure mechanisms [7] would be interesting studies towards the life-cycle prediction of such anomalous bio-tissues. Finally the variation of tissue properties by anatomical location, orientation, and layers of the UB motivates further studies on distributed order models incorporating nonlinearities. In this sense, mathematical and computational frameworks employing distributed-order viscoelastic models that learn distributions for fractional orders would be able to account for multi-fractal heterogeneous media and stochastic effects, leading to predictive zero-dimensional models with a reduced number of parameters.