METHODS FOR MODELLING CHANGES IN VISCOELASTICITY OF THE URINARY BLADDER BY ANATOMICAL LOCATION AND SWELLING By Laura Alison Nye A THESIS Submitted to Michigan State University i n partial fulfillment of the requirements f or the degree of Mech anical Engineering Master of Science 2020 ABSTRACT METHODS FOR MODELLING CHANGES IN VISCOELASTICITY OF THE URINARY BLADDER BY ANATOMICAL LOCATION AND SWELLING By Laura Alison Nye Urinary bladder dysfunction affects millions worldwide. It adds burden t o the healthcare system and individual patients with surgeries and long - term treatments such as daily catheterization. Patient specific modelling has been shown to red uce healthcare costs. Computational models are based on mechanical properties of tissues and can help in diagnoses or treatment plans. Although much work has been done on organs such as brain and heart, work on the urinary bladder is scarce. The bladder is a complex organ that exhibits time dependent behavior and so many factors must be consid ered when studying its mechanical properties. pathological causes of mechanical behavior, and modelling the viscoelasticity. We focus on two issues that will improve computational models of the urinary bladder wall. The first issue is to identify differences in mechanical behavior based on anatomical location and bath osmolarity. The second issue is to find an appropriate viscoelastic constitutive eq uation. Parameter characterization of viscoelastic models is especially challenging due to the time dependence of certain parameters. We explore the methods applied in literature and propose four possible models that would be appropriate for our experiment . The models reveal that for best results, we must normalize our data, choose an appropriate relaxation spectrum that has a unique solution, and consider nonlinear el asticity in addition to viscoelasticity. Preliminary results from these models suggest th at the lower body and trigone regions of the bladder have lower compliance than other regions. Our models also indicate a change in compliance based on bath osmolarity . In the future we will improve these results for definitive parameters that can be compa red statistically. We will do this through the implementation of nonlinear elastic viscoelasticity and triphasic theory. iii ACKNOWLEDGEMENTS First and foremost, I would like to thank my advisor Dr. Sara Roccabianca, who has known me since I was an undergrad learned much from her over the past few years and a m very thankful for all the time she has spent giving me guidance, reviewing my work, and being supportive overall. I would also like to thank others who made this project possible : m y thesis committee Dr. Tamara Reid Bush and Dr. Lik - Chuan Lee for their support , the creators of FEBio Gerard Ateshian and Steve Maas who answered many of my tedious questions, and Stacy H o llon in the Mechanical Engineering Department for helping me navigate many issues with patience and kindness. My lab members Mar issa Grobbel, Yuheng Wang, Dr. Sheng Chen, Tyler Tuttle, Dr. Mayank Sinha, and Eli Broemer supported me throughout this degree a s well with brainstorming sessions and much needed breaks in the computer lab. Lastly, I would like to thank my parents L inda Ny e and Paul Nye , who instilled the importance of education in me since a young age and without whom this would not have been poss ible. iv TABLE OF CONTENTS LIST OF TABLES ................................ ................................ ................................ ................................ ....... vi LIST OF FIGURES ................................ ................................ ................................ ................................ .... vii 1. RATIONALE AND MATHEMATICAL BACKGROUND ................................ ................................ . 1 1.1. Motivation : Urinary Bladder Dysfunction ................................ ................................ ...................... 1 1.2. Constitutive Model: Hyperelastic and Viscoelastic ................................ ................................ ........ 2 1.2.1. Continuum Mechanics ................................ ................................ ................................ ........... 3 1.2.2. Hyperelasticity ................................ ................................ ................................ ....................... 5 1.2.3. Donnan Equilibr ium Biphasic Model ................................ ................................ .................... 5 1.3. V iscoelasticity ................................ ................................ ................................ ................................ 6 1.3.1. The Weichert Model ................................ ................................ ................................ .............. 7 1.3.2. Linear Viscoelasticity ................................ ................................ ................................ .......... 10 1.3.3. Quasilinear Viscoelasticity (QLV) ................................ ................................ ...................... 11 1.3.4. The Reduced Relaxation Function ................................ ................................ ....................... 12 1.3.5. Advantages and Disadvantages to DRS and CRS ................................ ............................... 14 2. LITERATURE REVIEW ................................ ................................ ................................ ..................... 16 2.1. Urinary Bladder Physiolog y and Tissue Properties ................................ ................................ ...... 16 2.1.1. Submucosa ................................ ................................ ................................ ........................... 17 2.1.2. Detrusor Muscle ................................ ................................ ................................ ................... 18 2.1.3. Mechanical Properties vary by Anatomical Location and Orientation ................................ 18 2.2. Mechanical Tests ................................ ................................ ................................ .......................... 21 2.3. Experimental Challenges ................................ ................................ ................................ .............. 22 2.4. Challenges in Analytical Methods (Parameter Characterization) ................................ ................. 24 2.4.1. Intercorrela tion Algorithms ................................ ................................ ................................ . 25 2.4.2. Lack of Incremental Curve Fits ................................ ................................ ........................... 26 2.4.3. Inter - parameter Sensitivity as Mo tivation for Constant Relaxation Times ......................... 26 2.4.4. Method Selection ................................ ................................ ................................ ................. 28 3. MODELLING OF URINARY BLADDER TISSUE ................................ ................................ ........... 30 3.1. Introduction ................................ ................................ ................................ ................................ .. 30 3.2. Experimental Methods ................................ ................................ ................................ .................. 30 3.3. Previous Work ................................ ................................ ................................ .............................. 32 3.4. Optimizat ion Method ................................ ................................ ................................ .................... 34 3.5. Model 1: Quasilinear Viscoelasticity in FEBio, Single Stress - Relaxation ................................ ... 34 3.5.1. FEBio Model ................................ ................................ ................................ ....................... 34 3.5.2. Methods: Initial Results ................................ ................................ ................................ ....... 38 3.5.3. Methods: Guided Optimization for Quasilinear Viscoelasticity ................................ .......... 40 3.5.4. Results ................................ ................................ ................................ ................................ .. 41 3.5.5. Discussion ................................ ................................ ................................ ............................ 46 3.6. Model 2: Weichert Model with L ogarithmic Exponential Distribution ................................ ....... 46 3.6.1. Methods ................................ ................................ ................................ ................................ 46 3.6.2. Results ................................ ................................ ................................ ................................ ... 47 3.6.3. Discussion ................................ ................................ ................................ ............................ 49 3.7. Model 3: Reduced Relaxation Function with Log Normal Distribution (n=1) ............................ 50 v 3.7.1. Methods ................................ ................................ ................................ ............................... 50 3.7.2. Results ................................ ................................ ................................ ................................ .. 52 3.7.3. Discussion ................................ ................................ ................................ ............................ 53 3.8. Model 4: Reduced Relaxation Function with Log Normal Distribution (n=2) ............................ 54 3.8.1. Methods ................................ ................................ ................................ ............................... 54 3.8.2. Results ................................ ................................ ................................ ................................ .. 55 3.8.3. Discussion ................................ ................................ ................................ ............................ 57 3.9. Overall Conclusions ................................ ................................ ................................ ..................... 57 3.9.1. Sources of Er ror ................................ ................................ ................................ ................... 57 3.9.2. Future Work and Works in Progress ................................ ................................ .................... 58 BIBLIOGRAPHY ................................ ................................ ................................ ................................ ....... 61 vi LIST OF TABLES Table 1. Comparison of the Cauchy stress functions for Neo - Hookean elastic versus Neo - Hookea n hyperelastic materials as defined in the FEBio software (Maas et al., 2011) ................................ ............... 5 Table 2. Constants used in the Donnan equilibrium model of swelling used in FEBio (Maas et a l., 2011) 6 Table 3. P - values comparing statistical significance between peak stress es (left) and relaxed stresses (right) based on anatomical location ................................ ................................ ................................ ...................... 33 Table 4. P - values comparing statistical significance between peak stresses (left) and relaxed stresses (right) based on bath osmolarity ................................ ................................ ................................ ............................ 33 Table 5. Average fitted parameters by anatomic al location (top) and bath osmolarity (bottom) for the linear viscoelastic FEBio model (model 1) ................................ ................................ ................................ ........... 42 Table 6. Average RMSE values for model (1) ................................ ................................ ........................... 42 Table 7. Average fi tted parameters fitted to Eq. 43 and Eq. 44 by anatomical location (top) and swelling (bott om) for the logarithmic exponential CRS (model 2) ................................ ................................ ........... 47 Table 8. Average RMSE values for model (2) ................................ ................................ ........................... 48 Table 9. Parameters for Eq. 31 with n=1 fit to data for location t rials (top) and swelling trials (bottom) for model (3) ................................ ................................ ................................ ................................ ..................... 52 Table 10. Average RMSE values for model (3) ................................ ................................ ......................... 52 Table 11. Initial guess parameters from Nagatomi et al. (2008) ................................ ................................ 54 Table 12. Parameters for Eq. 31 with n=2 fit to data for location trials (top) and ba th osmolarity trials (bottom) for model (4) ................................ ................................ ................................ ................................ 55 Table 13. Average RMSE values for model (4) ................................ ................................ ......................... 55 vii LIST OF FIGURES Figure 1. Nonlinear viscoelasticity during incremental stress - relaxation experiment ................................ . 3 Figure 2. Reference (undeformed, material) configuration is transformed to current (deformed, spatial) configuration ................................ ................................ ................................ ................................ ................. 4 Figure 3. Voigt/Kelvin model (spring and dashpot in parallel ), the Maxwell model (spring and dashpot in series), and the Standard Linear Solid (SLS) (spring and Maxwell element in parallel) [Vincent, 2012] and their corresponding stress response to step strain [Fung, 1981] ................................ ................................ ... 7 Figure 4 . Weichert model schematic, comp osed of one spring element and Maxwell elements in parallel ................................ ................................ ................................ ................................ ................................ ...... 7 Figure 5. An example Prony Series on the logarithmic scale shows the decay of each component and contribution to initial value at ................................ ................................ ................................ ............ 10 Figure 6. Visual rep resentation of discrete spectrum of relaxation constants (top left) and how a relatively infinite spectrum of relation constants (top right) can be represented by a continuous spectrum (bottom) 13 Figure 7. An example of a DRS and CRS representing the same reduced relaxation function of time where [Shanbhag, 2019] ................................ ................................ ............................... 14 Figure 8. The geometry of the urinary bladder and layers of the urinary bladder wall [Roccabianca and Bush, 2016] ................................ ................................ ................................ ................................ ................. 17 Figure 9. The dissection plane on porcine UB used by Korossis et al. (left) and the anatomical locations used in their experiment (right) [Korossis et al., 2009] ................................ ................................ .............. 19 Figure 10. Relaxation spectrums of the ECM and SM [Nagatomi et al 2008] ................................ .......... 20 Figure 11. The overall stress relaxation function (left) determined by the equations on the right where and are constants determined by optimization [Nagatomi et al., 2008] ................................ ............ 20 Figure 12. Stress and strain of three common mechanical tests from left to right: stress - relaxation, creep, and dynamic ................................ ................................ ................................ ................................ ............... 22 Figure 13. Inter - dependence o f relaxation strength and relaxation time where , (corresponding to , in Eq. (21) ) change values when optimized to a new time constant ................................ ................................ ................................ ................................ ................................ .... 27 Figure 14. A UB sample mounted on the u niaxial machine, submerged in Krebs solution. ..................... 31 Figure 15. Five strains applied to the uniaxial machine (left) and the recorded stress relaxation behavior (right) ................................ ................................ ................................ ................................ .......................... 31 Figure 16. Samples cut for the location mechanical test, specifically dorsal (D), ve ntral (V), lateral (L), lower body (LB), and trigone (T) locations (left) and samples cut from the lateral location for the swelling test (right) ................................ ................................ ................................ ................................ .................... 32 Figure 17. Average peak stresses for first three strains 0.25, 0.5, and 1.0 for peak stresses (left) and relaxed stresses (right) for location trials ................................ ................................ ................................ ................. 32 Figure 18. Average peak stresses for first three strains 0.25, 0.5, and 1.0 for peak stresses (left) and relaxed stresses (right) for swelling trials ................................ ................................ ................................ ................ 33 viii Figure 19. Diagram of the optimization procedure ................................ ................................ .................... 34 Fig ure 20. Porcine UB mounted on the uniaxial machine (left) and an early FE model of the experiment showing stress in the x - direction (right) and the cross section where the forces were calculated as a dotted line ................................ ................................ ................................ ................................ .............................. 36 Figure 21. An isometric view of 1/8 th the origin al tissues with a 7x7x7 mesh (left) and the stress of the model immediately after the first strain where the color bar is stress [Pa] and surfaces A, B, and C are cross sections shown in Fig. 22 (right) ................................ ................................ ................................ ................ 37 Figure 22. FEBio model during a relaxation phase s howing 1/8 th of the original tissue at different cross sections (labelled in Fig. 21 above) ................................ ................................ ................................ ............ 37 Figure 23. Optimization by procedure in Fig. 23 (RMSE = 6.29) versus hand fitted results (RMSE = 4.35) and raw data on a normal time scale (left) and loga rithmic time scale (right) ................................ ............ 38 Figure 24. Relaxation times (pictured , ) are the location of the inflection point, on the logarithmic time scale, of the individual exponentials (boxes) and they decay to approxim ately zero within one logarithmic decade (bubbles) ................................ ................................ .... 39 Figure 25. The geometric mean between two time constants is the logarithmic hallway point where the decay from is complete if . ................................ .................... 40 Fi gure 26. Sections numbered in order that they are used for fitting parameters where section (1) is used for step (1), section (2) - (3) are used for step (2) where , and section (4) is used for step (3) ................................ ................................ ................................ ................................ ................................ .... 41 Figure 27. Average RRF (G(T)) pl otted for each anatomic al location dorsal (D), ventral (V), lateral (L), lower body (B), and trigone (T) for model (1) ................................ ................................ ............................ 43 Figure 28. Average RRF plotted for each osmolarity trial 100, 300, 600, and 900 mOsm/L and dry for model (1) ................................ ................................ ................................ ................................ ..................... 43 Figure 29. Average str ess relaxation curves p lotted for each anatomical location dorsal (D), ventral (V), lateral (L), lower body (B), and trigone (T) for model (1) ................................ ................................ .......... 44 Figure 30. Average stress relaxation curves plotted for each osmolarity trial 100, 300, 600, and 900 mOsm/L an d dry for model (1) ................................ ................................ ................................ ................... 45 Figure 31. The incremental stress - relaxation curves for location trials using parameters from Table 5 for model (2) ................................ ................................ ................................ ................................ ..................... 48 Figure 32. The incremental stress - relaxation curves for swelling trials using parameters from Table 5 for model (2 ) ................................ ................................ ................................ ................................ ..................... 48 Figure 33. Experim ental versus theoretical data for a good fit (bladder #2 900 mOsm/L) (left) and a bad fit (bladder #2 dry) (right) using model (2) ................................ ................................ ................................ ..... 49 Figure 34. Peak and relaxed stresses for each stress relaxation curve are used to normalize the cur ve by Eq. 45 [Sarver et al., 2003] ................................ ................................ ................................ ......................... 50 Figure 35. Normalized for all strains in the same experiment bladder #6 dorsal (left) and bladder #3, 100 mOsm/L (right) show no trends in increasing final stress that are seen in the nonnormalized data .... 51 Fig ure 36. Average relaxation behavior (RRF) by location based on param eters from Table 7 for model (3) ................................ ................................ ................................ ................................ ................................ 52 ix Figure 37. Average relaxation behavior (RRF) by bath osmolarity based on parameters from Table 9 for model (3) ................................ ................................ ................................ ................................ ..................... 53 Figure 38. Guided optimization for Log - proc edure outlined in Fig. 19 ................................ ................................ ................................ ...................... 54 Figure 39. Average RRF based on parameters from Table 12 for location trials for model (4) ............... 56 Figure 40. Average RRF based on parameters from Table 12 for swel ling trials for model (4) ............... 56 Figure 41. Viscoelastic solid with EFD neo - Hookean elastic component optimized in MATLAB (left) and manually fitted (right) ................................ ................................ ................................ ................................ . 58 Figure 42. FEBio triphasic model showing stress in the x - direction in different stages (A) loading t he fixed charge density of the Donnan E quilibrium neo - Hookean solid (B) loading the bath osmolarity (C) applying a strain to the specimen (D) the relaxation phase ................................ ................................ ....................... 59 Figure 43. The x - stress through the different stages of analysis of triphasic solid showin g higher peak stress for ramping bath osmolarity from 300 to 900 mOsm/L (left) and lower peak stress for ramping bath osmolarity from 300 to 100 mOsm/L (right) ................................ ................................ .............................. 60 1 1. RATIONALE AND MATHEMATICAL BACKGROUND 1.1 Motivation: Urinary Bladder Dysfunction Millions of peo ple suffer from urinary bladder dysfunction and th is condition has the potential to significantly decrease their quality of life and increase healthcare costs for them as individuals as well as for the US as a nation. Overactive bladder (OAB), undera ctive bladder (UAB), and neurogenic bladder are a few of the conditions requiring patients to accrue daily costs due to, for example, catherization or incontinence pads, and ultimately potentially major surgery. For example, one of the long - term solutions that a re used to treat urinary incontinence is augmentation surgery. This procedure can potentially significantly increase bladder functionality and, ultimately, the quality of life of patients affected by incontinence . H owever, there is a somewhat high risk ass ociated with it , in fact approximately 5% of patients undergoing this kind of surgery have repaired bladders rupture due to material incompatibility [Korossis et al., 2009]. Furthermore , Flack et al. (2015) pointed out that while surgical inte rvention may reduce costs long - term for an entire healthcare system, it is not usually financially feasible for individual patients. It was estimated in 2007 that the total cost of urinary incontinence in the United States was $65.9 billion for adults 25 y ears and olde r. There is potential, however, for reducing this cost with patient specific modelling. This is shown, for example, by a study that used portable bladder ultrasounds to reduce catheter costs by 46% based on a nt [ Polliack et al., 2005 ]. This suggests that personalized medicine has great potential to decrease healthcare costs and increase the quality of care for each patient. To this end, in the last twenty years, other fields like cardiology and neurology have used mathemat ical modeling to generate patient - specific models that could significantly help clinicians with diagnosis and treatment choices. A widely used modeling technique, that has showed some level of success in other organs in the body such as the he art, brain, a nd skin, is finite element modelling (FEM) based on realistic mechanical description of tissues [Roccabianca and Bush, 2016]. The potential of these models is that they consider 2 how mechanical behavior of tissues changes due to age, sex, disea se, injury, m icrostructure , and organ geometry, all of which have been identified as important factors in whole organ function. Despite the enormous amount of literature on biomechanical modelling of soft biological tissues, work focusing on the urinary b ladder (UB) h as been largely overlooked. More accurate models are needed to predict the feasibility of tissue engineered scaffolds, preoperative planning for surgery, and preventive care. Theoretical modelling of urinary bladder tissue based on experimenta l data has a great potential to help diagnose disease, improve tissue engineered materials [Dahms et al., 1998], analyze interactions with medical devices [Natali et al., 2015], and ultimately contribute to the reduction of healthcare costs associated wit h urinary dys function. 1.2 Constitutive Model: Hyperelastic and Viscoelastic It is well established that UB tissue has viscoelastic and hyperelastic behavior. Both v iscoelasticity and hyperelasticity ha ve been modelled in many papers on the UB [Fry and Wagg, 1999] [Nag atomi et al., 2004] [Nagatomi et al., 2008] [ Van Mastrigt and Nagtegaal, 1981] [Thiruchelvam et al., 2003] [Natali et al, 2015]. Hyperelastic models are suitable to describe experimental data collected under pseudo - static conditions and are incapable of de scrib ing stress - relaxation experiments on their own. Similarly , linear viscoelastic models are only suitable to describe strains less than 0.1 [Vincent , 2012 ] and do not model stress - relaxation of the UB sufficiently. This can be observed i n stress - relaxat ion experiments involving multiple strain increases. In a stress - relaxation test of UB tissue, a fter the initial strain , a sharp decline in stress is observed, followed by a time dependent stress reduction until the sample becomes fully rel axed ( Figure 1 ). At equilibrium , the elastic component of the tissue is the only contributor to stress, according to quasilinear viscoelasticity ( QLV ) theory [Fung, 1981]. If the elastic component behav ed in a purely linear way , a second strain would produ ce a stress at e quilibrium that is linearly proportional to the strain increase. However, t his is not the 3 case in our experimental data ( Figure 1 ). Therefore, viscoelasticity must be paired with a nonlinear elastic component rather than a linear elastic one. Figure 1. Nonlinear viscoelasticity during incremental stress - relaxation experiment In the remaining sections of Chapter 1 I will introduce the mathematical framework that has been used in the past to study both elastic and viscoelastic be haviors of soft biological tissues. First, I will briefly discuss h yper e lastic models and biphasic models , then I will focus in greater detail on discussing v iscoelastic models, as it was the focus of parameter characterization difficulties. 1.2.1. Continuum Mechanics In the conti nuum mechanics framework, one can use the deformation gradient tensor, , to represent the displacement from the reference configuration (undeformed) to the current configuration (deformed) . The location of each particle in the reference and current configurations are described the position vectors (material de scription) and (spatial description), respectively ( Figure 2 ). The displacement vector u is the difference between the two: . 4 Figure 2. Reference (undeformed, material) configuration is transformed to current (deformed, spatial) configura tion The deformation tensor relates the material and spatial coordinates , and is defined as : or . (1) From the deformation tensor, one can evaluate the r ight Cauchy - Green d eformation t ensor ( , material de scription ) and the l eft Cauchy - Green d eformation t ensor ( , spatial description ) , as : , (2) . (3) Furthermore, one can identify two commonly used t wo strain tensors , such as the Green - Lagrange t ensor ( , material description ) and the Alamansi s train t ensor ( , spatial description ) defined as: , (4) (5) w here I is the identity matrix. The deformation metrics described above are used in many constitutive equations, such as neo - Hookean material description, Donnan equilibrium model, and quasi - linear viscoelasticity. 5 1.2.2. Hyperelasticity Soft biological tissues have nonlinear relationship s between stress and deformation and often underg o large deformations ; for these reasons they re quire more sophisticated constitutive equations than linear elasticity to describe their behavior . There is no consensus on the best hyperelastic model to describe the pseudo - static behavior of the urinary bladder ; in our study we first employed the neo - Ho okean hyperelasticity , in order to reduce the number of variables. Table 1 compares this model wi th where and and are the Lam é parameters defined as: (6) (7) where and Table 1. Comparison of the Cauchy stress functions for Neo - Hookean elastic versus Neo - Hookean hyperelastic materials as defined in the FEBio software (Maas et al. , 2011 ) 1.2.3. Donna n E quilibrium B iphasic M odel There are multiple approaches to modelling swelling of tissues in an aqueous solution. In the FE software FEBio, for example, users can opt to use biphasic - solute, biphasic, triphasic or multiphasic theory, and combine these wi th different types of solutio ns, porosity, and material models. The biphasic model describes a mixture of a porous permeable solid and interstitial fluid (Maas et al. , 2011 ). To model a triphasic material with the Donnan equilibrium model of swelling, one assumes that a negative charg e 6 gradient between the solution and the solid prompts movement of fluid into or out of the porous solid, resulting in swelling or shrinking of the tissue. Th e inputs of the model are the bath osmolarity , as well as elastic prop erties, the volume fraction, and the fixed charge density of the porous media . The effective stress of the mixture is a combination of elastic stress and fluid pressure : , (8) where the fluid pressure is defined as : (9) The constants introduced in Eq. 9 are dependent on the properties of the solid, fluid, and surroundings, and are specified in Table 2 . Table 2. Constants used in the Donnan equilibrium model of swelling used in FEBio (Maas et al., 2011) 1.3. Viscoelasticity The focus of this thesis is to identify an accurate model to d escribe the viscoelastic properties of the UB wall . Mathematical models that describ e viscoelasticity derive from mechanical spring - dashpot models. The three basic linear viscoelastic models found in literature are the Voigt /Kelvin model, the Maxwell model , and the Standard Linear Solid (SLS), pictured in Figure 3 . 7 Figure 3. Voigt/Kelvin model (spring and dashpot in parallel), the Maxwell model (spring and dashpot in series), and the Standard Linear Solid (SLS) (spring and Maxwell element in parallel) [Vincent, 2012] and their corresponding stress response to step st rain [Fung, 1981 ] 1.3.1. The Weichert Model Complex arrangements of elements of springs and dashpots in series or parallel (known as lumped parameter models) have been explored by Tschoegl (1989) . However , literature suggests that the Weichert m odel ( also called generalized Maxwell model) is acceptable for material characterization of viscoelastic soft tissues [Vincent , 2010 ]. This model is composed of a single spring element and Maxwell elements in par allel ( Figure 4 ). Figure 4. Weichert model schematic, composed of one spring element and Maxwell elements in parallel 8 T he Maxwell model is developed under the following assumptions : (1 0 ) ( 11 ) Where , , are from Figure 3 . The strain and strain rate generate stresses in the spring and in the dashpot elements, respectively , following the laws: ( 12 ) ( 13 ) where k is stiffn ess and is viscosity . Combining equations (1 0 ) ( 13 ) one can calculate a relationship between strain rate , and stress and stress rate, as: ( 14 ) Finally, d efining the time c onstant as and rearranging equation ( 14 ), we can ultimately derive the differential equation for stress : ( 15 ) When this model is applied to a stress relaxation test, strain rate is assumed to be 0 after applying the stress. In serting and integrating, the final time dependent stress can be written as: ( 1 6) From Eq. (16) it can be seen that for the Maxwell model when , the stress goes to zero. However, soft tissues do not relax to a st ress equal to zero when , but they relax to a value of stress dependent on the applied strain. This is why the Maxwell model is combined in parallel with a simple spring in the SLS model . It is easiest to derive the SLS stress - strain relationship from the Laplace transform of Eq uations (12) and (16) , where . The stress of the Maxwell and spring elements, respectively, are then: 9 ( 1 7) ( 1 8) Using the assumption that , t he stress in the Laplace space can be written as : ( 1 9) w here is the time depende - relaxation test, strain is zero until time and can be represented by the Heaviside function. For , we have that : ( 20 ) Inserting this expression for strain into Eq. (19) and taking the inverse Laplace yields the final stress - strain relationship for the Weichert Model , with Maxwell elements , as: ( 21 ) In Eq. (21) and are referred to as the relaxation strengths and relaxation times, respectively ; is representative of the modulus of elasticity of the elastic ( spring) element; and is the relaxat ion function that defines t he discrete spectra in pairs . The summation in Eq. (21) is also referred to as the Prony Series . A visual representation of how each component of the series contributes to the overall behavior is seen in Figure 5 . 10 Figure 5. An example Prony Series on the logarithmic scale shows the decay of each component and contribution to initial value at 1.3.2. Linear Viscoelasticity The Weichert Model is a good starting point for understanding the derivation of viscoe lastic models but is usually applied only to 1D uniaxial tension tests. Additionally, derivation of the stress - strain relationship from the Laplace domain to time domain requires knowing the strain history, t herefore, a general model of viscoelasticity app licable to stress and strain tensors with any strain history is desired. A one - dimensional case of a simple bar elongated by and experiencing stress is considered. All the lengthening history due to is considered up until time t . For a small time - interval , the elongation increment is and causes a stress . The elongation and st ress increments are related to each other via the proportionality constant which is a funct ion of the time interval ( ) as : ( 22 ) It follows that we can estimate the entire loading history by integrating Eq. ( 22 ) and obtaining : 11 ( 23 ) where is the relaxation function. If is assumed to be a tensor of material characteristics, then the stress - strain relationship in Eq. ( 23 ) can be written as a constitu tive equation, with and as function of space and time, and bein g the tensorial relaxation function : ( 24 ) If motion begins at , and for , then Eq. (24) becomes ( 25 ) A viscoelastic solid that behaves following the law described in Eq. ( 25 ) is s aid to be linear viscoelastic. However, tissues have nonlinear stress - strain relationships and so Fung (1981) developed a different viscoelastic function . In continuum mechanics, the QLV model has broader applications because it can be used with most exper imental data sets such as stress relaxation, creep, and hysteresis [Fung, 1981]. QLV theory has been applied in finite element softw are, including the viscoelastic model found in FEBio [Maas et al., 2011]. 1.3.3. Quasilinear Viscoelasticity (QLV) Consider a cylindrical specimen undergoing an infinitesimal change in stretch where the stress response is assumed to be a function of time t and stretch . The stress response denoted as , is also known as the relaxation function and is defined as : ( 26 ) where is the reduced relaxation function (dimensionless) and is the elastic response (dimension of a stress). The stress response for an instant of time where can be evalu ated as : 12 ( 27 ) It follows, that the stress a t any time t , which is depending on the entire stress history, can be evaluated from Eq. ( 26 ) as : ( 28 ) Again, the strain is assumed to be applied at , as well for , which allows us to evaluate : ( 29 ) which can be generalized in a 3 - dimensional c onstitutive equation as : ( 30 ) where is the Kirchhoff stress tensor, is the reduced relaxation function tensor. , or that the function is nondimensional. This function resembles Eq. ( 25 ) except for the elastic stress can be derived from any hyperelastic or ela stic material. 1.3.4. The Reduced Relaxation Function At we assume that . This is accomplished with the normalization of the relaxation spectrum by: (31) The purpose of normalizing the relax ation spectrum is so that the relaxation parameters can be studied separate from elastic parameters. The condition implies that . This holds true for 13 normalized data, but it should be noted that in the Weichert model, do not su m to 1 because they are derived a generalized reduced relaxation function compatible with continuous spectrums : (32) The spectrum considers the variable to be continuous from 0 to . The relative amplitude (relaxation strengths in a discrete spectrum) is a function of (see Figure 6 ). The generalized QLV model ( E q. 30 ) can use either a d iscrete relaxation spectrum (DRS) or a continuous relaxation spectrum (CRS) where a DRS is a sum of discrete spectra pairs and the CRS takes the form of an integrable function of Figure 6. Visual representation of discrete spectrum of relaxat ion constants (top left) and how a relatively infinite spectrum of relation constants (top right) can be represented by a continuous spectrum (b o ttom) 14 1.3.5. Advantages and Disadvantages to DRS and CRS The same reduced relaxation function can often be repr esented equally satisfactorily by both CRS and DRS as seen in Figure 7 . However, they each come with benefits and downsides. Choosing which one is best can depend on the application. Figure 7 . An exam ple of a DRS and CRS representing the same reduced rel axation function of time where [Shanbhag, 2019] DRS are often used because of their computational efficiency [Park and Kim , 2001 ]. The equation is straightforward: it is the sum of exponential functions, referr ed to as the Prony Series in many papers. Unfortunately, the Prony series is an ill - posed problem, meaning it does not have a unique solution [ Baumgaertel and Winter , 1989 ] [McDougall et al. , 2 014 ] . In other words, a single relaxation curve can be represen ted by multiple combinations of relaxation spectra [Fung , 2013 ]. Additionally, DRS can include as many exponentials as desired for a well - fit curve, adding difficulty to the model selection bec ause the user must define n , the number of viscoelastic pairs . For example, Liu, Z., and K. Yeung (2008) used an RRF composed of three sets of exponentials , resulting in eight total parameters to be optimized to describe their experimental data. On the one hand, while i ncreasing the order of the Prony Series will decr ease the apparent fitting error , it will also increas e - Baumgaertel and Winter , 1989 ]. On the other hand, a lower order Prony Series will decrease the uncertainty but might increase the error between theory and experiment s . 15 A CRS comba ts this by reducing the number of variables. A CRS can reveal useful information about the relaxation behavior of constituents, as seen in the work of Nagato mi et al. (2008). Recent works have emphasized the need for rheologists to include CRS in their res ults [ McDougall et al. , 2104 ] and others have opted to only use them in their analyses of soft tissues [Nagatomi et al. , 2004 ] [ Nagatomi et al. , 2 008]. CRS can require as few as 3 variables to define the RRF, depending on the spectrum chosen. Although it reduces uncertainty, CRS are more difficult to implement because they involve integration. As a result of its computational expense, CRS are left out of vis coelastic models available in popular FE software. 16 2. LITERATURE REVIEW 2.1 Urinary Bladder Physiolo gy and Tissue Properties The UB is a multi - layered storage organ that undergoes large deformations multiple times a day, throughout the filling and voiding cycles. The urinary bladder wall (UBW) has complex mechanical behavior, that can be considered hyper elastic, in certain loadi ng conditions, and viscoelastic (time dependent), in others. Sudden increases in volume due to fluid entering the bladder rapidly increase the pressure followed by a phase of viscoelastic relaxation until reaching a steady state ag ain [Nagatomi et al., 200 4]. The mechanical properties of this organ are crucial to its function, indeed a change in the pressure - volume relationship of the wall due to disease can significantly alter the UB function. Furthermore, disease and injury can ca use remodeling of the UBW compliance, contractility, and elasticity. A study on fetal bladder outlet obstruction (BOO), for example, found that obstructed bladders exhibited increased compliance and decreased elasticity and viscoelasticity [ Thiruchelvam et al., 2003 ]. Likewise, contractile abnormalities can cause outflow obstruction [Fry and Wagg, 1999]. mechanical properties ar e tied to its physiology. Obstructed bladders exhibit collagen deposition and ultimately to increased stiffness, compliance [Fry and Wagg, 1999]. Furtherm ore, the bladders collected from rats subjected to a spinal cord injury were found to have increased elastin density, which resulted in reduced stress relaxation behavior (retaining a higher stress for longer) and increased bladder compliance [Nagatomi et al., 2004]. Therefore, it is important to consider the physiology of the UB and contribution of its components to its function. The urinary bladder wall is a complex structure that can be divided in four distinct layers, each with their own unique properti es: mucosa, submucosa, detrusor muscle, and adventitia ( Figure 7 ). The innermost layer of the UBW, called the mucosa, is a thin layer of epithelial cells , which gathers in folds when the 17 bladder is empty. Umbrella cells and a layer of glycosaminoglycans pr otect the rest of the UBW by blocking diffusion of urine from the lumen. The outermost layer, the adventitia, is a thin layer composed of connective tissues. The two thickest layers, the submucosa and detrusor muscle, bear most of the load during deformati on and have been the focus of mechanical testing on the UBW. Figure 8 . The geometry of the urinary bladder and layers of the urinary bladder wall [Roccabianca and Bush, 2016] 2.1.1. Submucosa The submucosa is capable of large def ormations, crucial to the filling phase, and consists mostly of extracellular matrix, that contai ns a network of collagen and elastin fibers. Elastin fibers provide compliance whereas collagen fibers provide strength. It is generally accepted that the inte ractions between fibers and the uncoiling and rearranging of fiber networks contribute to the mec hanical behavior, both elastic and viscoelastic, of soft tissues [Shen et al., 2011] [Natali et al., 2015]. Overall , the ECM has observable viscoelastic proper ties , as seen in experiments performed on isolated ECM [Nagatomi et al., 2008]. 18 2.1.2. Detrusor Muscle The detrusor muscle mechanical function is to relax and distend during filling and perform active contractions during emptying. Smooth muscle cells cont then drives voiding. Contractility of the bla dder smooth muscle cells (SMCs) is dependent on force generation within the muscle cells as well as viscous interactions with the ECM they reside in [Fry and W agg, 1999]. Determining whether abnormalities in contractility are due to the muscle cells themse lves or to changes in composition and structure of the ECM motivated Fry and Wagg to study isolated detrusor muscle in obstructed and healthy bladders [Fry and Wagg, 1999]. Obstructed bladders exhibit higher deposition of extracellular material compared to healthy bladders and offer a way to determine whether changes in the ECM are responsible for reduced compliance. They found no statistical difference in visco elastic properties of the detrusor muscle in the passive state in healthy versus obstructed bladd ers. However, in the active state, viscoelastic properties were more apparent in the obstructed bladders than the healthy bladders. They speculated that this c because of deposition of ECM material. 2.1.3. Mechanical Properties Vary by Anatomical Location and Orientation Multiple constituents and layers contribute to the macroscopically observed mechanical be havior in the UBW. It is difficult experimentally to separate the layers of the UBW, so often its mechanical properties are determined from bulk behavior of strips of the UBW containing all four layers. In this state, experimentalists can still observe ani sotropy due to fiber orientation. Uniaxial tests over the years have shown that mechanical pro perties differ between strips in the longitudinal and transverse directions as well as between anatomical locations. Korossis et al. (2009) tested porcine urinary bladder strips in longitudinal and circumferential directions in 5 locations of the bladder: dorsal, trigone, lateral, ventral, and lower body ( Figure 8 ). Uniaxial tests of the specimens were accompanied by histological characterization . They found that s trips oriented longitudinally had higher stiffnesses than those oriented circumferentially. Th is is supported by their histological results which shows that the elastin fibers are oriented primarily in the 19 transverse direction. There were insignificant dif ferences in anisotropy in the transverse direction by anatomical location, finding that is sup ported by the observed uniform expansion circumferentially during the filling phase. Regional differences were significant as well. The dorsal and ventral region s had the highest compliance and least directional anisotropy. The lateral, lower body, and tr igone regions had the highest tensile strengths, and the trigone region had the lowest failure strains along with the highest collagen phase slope, making it the least distensible region. In general, the basal regions (lower body and trigone) were less com pliant than the apex regions (dorsal and ventral) and these properties were attributed to their collagen and elastin networks. From their results, it can be concl uded that material anisotropy, elasticity, and viscoelasticity will vary based on harvest loca tion. Figure 9 . The dissection plane on porcine UB used by Korossis et al. (left) and the anatomical locations used in their experiment (right) [Korossis et al. , 2009] Other researchers focused on isolating constituents to observe their individual mechan ical properties, as opposed to study the tissue intact. As mentioned previously, Fry and Wagg (1999) studied the detrusor muscle independent of the rest of the UB W by dissecting detrusor muscle strips and observed significant viscoelastic properties. Nagat omi et al. (2008) isolated the ECM through decellularization in order to determine a reduced relaxation function (RRF) (derived in Section 1.3.4 ). for the ECM ind ependent of 20 the smooth muscles (SMs) ( Figure 9 ). They used a dual Gaussian spectrum ( Eq. 31 ) where , , , are the amplitude, mean, and variance, respectively. Additionally, is an integer that rep resents the total number of relaxation mechanisms are considered, and was chosen to account the mechanical behavior of both collagen and elastin (see dotted lines in Figure 9 ) . Results from a previous study on SMCs were used for a constant value conti nuous spectrum ( Eq. 32 ), where and are the start and end of the spectrum respectively and is the index representing the amount of overall relaxation and is dimensionless. ( 31 ) ( 32 ) Figure 10 . Relaxation spectrums of the ECM and SM [Nagatomi et al ., 2008] (33) ( 34 ) ( 35 ) Figure 1 1 . The overall stress relaxation function (left ) determined by the equations on the right where and are constants determined by optimization [Nagatomi et al., 2008] The overall UBW reduced relaxation function ( Fig . 10 , Eq. 34 ) was then modelled as the summation of each reduced relaxa tion function and determined by Eq. (33) multiplied by their respective recruitment functions and ( Eq. 35 ). (The physical meaning and mathematical 21 description of the RRF and relaxation spectrums was deta iled in Section 1.3.4. ) These two studies demonstrate that mechanical propert ies vary between constituents. While the exact contribution of individual components to viscoelasticity and hyperelasticity remains understudied, constitutive models are useful fo r predicting response to pathologies that influence geometry, composition, an d therefore the capabilities of the bladder to void urine. Regardless of the type of constitutive equation chosen to describe the mechanical behavior, for the description of the t issue to be accurate the constitutive parameters must be informed by experime ntal data obtained from mechanical tests. This is usually achieved by an optimization process, via error minimization. In the following, we will first describe the mechanical test s that have been previously used to quantify the behavior of the UBW, and sec ond the methods that have been used to fit models to experimental data. Here we focus on viscoelastic behavior. 2.2. Mechanical Tests To study viscoelastic behavior of soft biolog ical tissues, three common mechanical tests are employed: stress - relaxation, creep, and dynamic tests [Roylance, 2001] seen in Figure 11 . Stress - relaxation tests are usually performed on a uniaxial or biaxial te nsile machine and consist of applying an inst antaneous strain to the sample followed by a relaxation time at constant strain. The stress - time curve will show a decrease - ss - relaxation tests, but aim to maintain a co nstant stress over time, rather than a constant strain, this results in an increase in strain over time. Finally, dynamic loading tests are performed by loading a tissue surface with an oscillatory load and quan tifying the phase shift between applied stres s and recorded strain. Uniaxial stress - relaxation experiments were performed in this work. 22 Figure 1 2 . Stress and strain of three common mechanical tests from left to right: stress - relaxation, creep, and dynam ic 2.3. Experimental Challenges Replicating and comparing experimental results on the UB has presented many obstacles over the years. Results from studies often contradict each other. For example, two studies comparing the maximum stress of young versus ol der bladders found opposite relationships between age and stiffness. Martins et al. (2011) tested female cadavers for their maximum stress an d found that the group of females age 50 years and older had lower stiffness compared to the group of females below 50 years of age; they concluded that bladder stiffness decreased with age. Chantereau et al. (2014), however, in a similar study on female c adavers found that bladders harvested from older individuals, with an average age of 75, had a higher stiffness com pared to that of younger individual, with an average age of 29 . Another set of studies found conflicting results about the anisotropy of the UBW. While one study found a larger stiffness in the longitudinal (apex to base) direction [Korossis et al., 2009], the other study found higher stiffness in circumferential (transverse) direction [Zanetti et al., 2012]. Details of these studies results an d other UB experiments can be found in literature [Roccabianca and Bush, 2015]. These discrepancies could be due to a lack of standardized experimental protocol. 23 Studies often vary in strain rate, temperature, reference frame, and maximum strain applied ( high strain >200% or lower strain <100%). Many tests are conducted at room temperature and temperature is known to affect viscoelasticity [Roylance , 2001]. Papers also fail to specify important factors about their samples, such as sex [Dahms et al., 1998]. Lastly , experimentalists are challenged with creating ex vivo environments that mimic as much as possible in vivo conditions. For example, hydration seems to greatly affect the elastic and viscoelastic behavior of the UBW , however, specimens are often tested at different levels of hydration. Many papers do not mention submerging their specimens in any solution [Chant e reau et al., 2014], or they store samples in solution prior to testing but not during [Natali et al., 2015] [Dahms et al., 1998]. Others perform tests in modified Krebs solution [Nagatomi et al., 2004] [Nagatomi et al., 2008] [Van Mastrigt and Nagtegaal, 1 981], Tyrode solution [ Thiruchelvam et al., 2003 ] or saline solution [Korossis et al., 2009]. Furthermore, it has been suggested that viscoelasticity of tissues in vivo can be attributed, at least in part, to viscous interactions between cells and interst i tial fluid [Natali et al., 2015]. The large variation between experimental approaches makes it crucial to question the validity of comparisons between studies. For now, the findings indicate the complexity of modelling the UB based on a multitude of fact o rs: geometry, regional differences, fiber orientation, constituents, and experimental approach. Multiphysics models might someday model multiple phenomena at once during the filling and voiding of the UB. Until then, it is important to isolate individual f actors and study their effects on overall UB function to aid in the creation of computational models. In this study we focused on two factors of interest: the differences in viscoelastic behavior by anatomical region and by bath osmolarity during a swelli n g test. Work by Korossis et al. (2009) highlighted the inhomogeneity among anatomical region and direction of the bladder. From uniaxial quasi - static mechanical tests on porcine bladders, they found that the basal region of the UB around the neck of the u r ethra had higher stiffness than apical regions. To investigate viscoelastic behavior based on anatomical region and bath osmolarity, we must discuss the different viscoelastic models available. As discussed in Section 1.3 , a discrete relaxation spectrum ( D RS) or continuous relaxation spectrum (CRS) can be used. Choosing the best one can be difficult. 24 2.4. Challenges in Analytical Methods (Parameter Characterization) Parameter characterization of viscoelastic materials can be split into two parts. The elasti portion of the constitutive equation can be determined by the equilibrium stress at the end of a stress relaxation curve. This correlates with the sin gle spring element in the Weichert model ( Figure 4 ). The relaxation spectra are more difficu lt to approximate. optimization algorithms. The solutio n is not very straightforward due a few complications: 1) CRS and DRS are ill - posed problems me aning they do not have a unique solution. Optimization of spectra might arrive at false minima before finding the best fit spectra [Doehring et al., 2004 ]. 2) CRS and DRS are very sensitive to noise in the data [ McDougall et al. , 2014 ] 3) DRS optimized without c onstraints can give negative parameters that have no physiological meaning [ McDougall et al. , 2014 ] 4) QLV does not account for strain - rate dependence in large deforma tions [Fung, 1981 ] Of the first attempts to fit exact relaxation spectra, Schapery (1962) fi t parameters using a simple collocation but failed to prevent negative coefficients. To tackle negative parameters, later methods attempted interactive ad justment of parameters, recursive algorithms, regularization, and quadratic programming [Park and Kim, 2001]. Simpler approaches have involved pre - smoothing noisy d ata with a power series [Park and Kim, 2001]. Detailed discussion of these techniques in ter ms of computational efficiency, mathematical background, and implementation are readily available in literature [Park and Schapery , 1999 ] [Park and Kim , 2001 ] [Shanb h ag , 2019 ] [McDougall et al. , 2014 ]. And so, while issues (1), (2), and (3) are solved by t hese methods, they are relatively inaccessible. While many have published papers of their work fitting CRS, very few have shared their codes online , as pointed out by Shanbhag and McDougall et al. in their works [ Shanbhag, 2019 ] [ McDougall et al., 2014 ]. E nd - users of these codes are likely to be unfamiliar 25 with mathematical concepts presented in literature making it difficult to reproduce code independently. Rec ently, computational physicists have created open - [Shanbhag , 2019 ]. To address strain - rate dependence, issue (4), QLV has been altered on a few occasions in literature [ Van Mastrigt and Nagtegaal , 1981 ]. Ab ramowitch and Woo (2004) developed a piecewise constitutive equation to simultaneously fit the ra mping and relaxation phases. It was shown to improve the fitting of experimental data sets in comparison with QLV theory alone. Furthermore, more complications arise depending on the form of stress - strain data and comparisons one wishes to make between the data sets. Finally , these methods have not been applied to 3D constitutive equations , to our knowledge. The major concerns are then: 5) CRS cannot be used to sim ulate tissues and organs in FE 6) Algorithms available online do not fit data from multi - step strain s and each stress - strain curve must be analyzed independently QLV requires normalizing stress data to the peak stress 7) Inter - parameter sensitivity (see Figure 13 ) makes it difficult to compare fitted data among samples differing in pathology Our research ai ms to fit data to 5 incremental strains resulting in a series of stress - relaxation curves ( Figure 1 ). Additionally , we hypothesize based on work by Korossis e t al. (2009) that there may be differences in viscoelastic parameters based on anatomical location s. Therefore, we must be able to compare parameters between the data sets. Implementing these models in FEBio is the final goal. 2.4.1. Intercorrelation Algor ithms A solution to the issue (5) has been addressed by intercorrelation algorithms that solve for the CRS and create DRS based on the CRS. Intercorrelation was successfully used by Labus et al. (2016) to model brain tissue by using peaks in the CRS to def ine the best discrete relaxation times. Also, Mun and Goangseup (2010) used an intercorrelation me thod to fit 16 pairs of relaxation constants to CRS obtained from dynamic 26 intercorrelation technique to solve for the DRS based on CRS. 2.4.2. Lack of Incremental Curve Fit s It is standard practice to fit elastic parameters to raw stress data, normalize the stress data, then fit viscoelastic parameters to the normalized stress. This method separates elastic and viscoelastic parameters, which are meant to be dimensionless [Sa rver et al. , 2003 ]. This works well for single stress relaxation curves but was rarely applied to incremental stress - relaxation curves for experiments with se quential strain The stress state of the tissue is dependent on all previously applied loads. Many papers avoid fitting incre mental tests by only ramp ing their samples once. Liu and Yeung (2008) fit their parameters to normalized stress curves of s trains of 5%, 10%, and 15% separately. The highest strain 15% had the lowest relaxed normalized stress. This is contradictory to our raw data where incrementally larger strains were applied to samples and larger strains produced larger r elaxed stresses. Sa rver et al. (2003) comment on this relationship between strain and relaxation; they note that common normalization er et al. (2003) de rive an expression for stress relaxation based on incremental strains using a non - strain dependent normalization method, solving issue (6). Unfortunately, this expression differs from constitutive equations used in FE software such as FE Bio, so parameters optimized will not transfer to 3D simulations. 2.4.3. Inter - parameter Sensitivity as Motivation for Constant Relaxation Times For situations where only discrete relaxation spectra can be used, such as in FE simulations, caution is advised for parameter co mparisons of bladders with different pathology. One may fit relaxation strengths and relaxation times to separate sets of data, but it is not wise to draw correlations between the relaxation strengths if the relaxation times are dif ferent . The relaxation s trengths are coefficients to the exponential, 27 which is dependent on the relaxation time . To demonstrate this, a Prony Series ( Eq. 21 ) of 3 exponentials was fit to stress relaxation data from a single strain from our experiments usin g the time constants , and a range of from to . The resulting values of change every time is increased ( Figure 13 ). The entire Prony Series is a summation of all relaxation pairs so any change to will require a chan ge to and as well ( Figure 13 ). Figure 13. Inter - dependence of relaxation strength and relaxation time where , (corresponding to , in Eq. (21) ) change values when opti mized to a new time constant Additionally, the DRS represented by a Prony series has sensitivity to different parameters at different location along the curve (see Figure 5 ) . Simple least squares optimization is inefficient because they do not accou nt for parameter relationships. This is likely why many choose to set their relaxation times as constants and optimize the relaxation strengths. [Goh et al. , 2004 ] fi t between two to five relaxation strengths to nonlinear elastic polynomial and Ogden model s. Abramowitch and Woo (2004) fit their piecewise model using three set time constants using the Levenberg - Marquardt approach to minimize their function. This is the best way to fit time constants with certainty in results if using a least - squares minimiza tion approach. The downside to this method is that the investigator must determine the best time constants and best number of time constants heuristically, introducin g human bias. 28 More methodical approaches avoid user bias by using machine learning instea d of setting the relaxation times as constants. The potential usefulness of machine learning is that it can solve ill - posed problems like the Prony Series without use r - bias. Genetic algorithms are a subset of machine learning techniques inspired from the D arwinian principle of survival of the fittest [ Kohandel , 2008 ] . The algorithm begins with a n possible solutions and allows them to progress to the next evolution phase if they are best fit compared to their competitors. Kohandel (2008) use d a genetic algorithm to fit both relaxation strengths and relaxation times to the alternative QLV constituti ve equation developed by Abramowitch and Woo (2004) . A genetic algorithm is available in MATLAB that considers relationships between parameters. Th is can be used to create constraints between variables using inequalities. 2.4.4. Method Selection In conclus ion, there exists algorithms developed over decades to fit exact relaxation spectra for DRS and CRS alike. These are incredibly useful for studying the relaxation behavior of soft tissues through a single step strain test. These algorithms are likely to ai d in identifying underlying mechanisms of relaxation in tissues in uniaxial stress - relaxation tests. If the only desire is to simulate bladder beha vior in FE software such as for patient specific modelling from cystography, CRS conversion to DRS is a suita ble option. In the case where relaxation behavior is to be compared between data sets rather than averaged from all tests, it is advisable to set r elaxation times for the DRS to avoid fitting parameters that are not actually comparable. A methodical approa ch for fitting relaxation strengths to set relaxation times will speed up analysis since sensitivity of the relaxation constants along the curve va ry. Fitting a continuous relaxation spectrum is also useful for comparison of these data sets but should not be converted to DRS because those will likely yield different relaxation times, leaving the parameters incomparable. For complex strain histories, researchers might need to replicate models such as those by Abramowitch and Woo (2004) and Sarver et al. (200 3) or use FEBio recognizing that time dependent and elastic components are inseparable for analysis. Parameters from the former method are not tran sferable to the latter method. 29 If FE simulation of complex strain histories are desired for comparison of dat a sets, as in our case, the user must set their relaxation times if they hope to find statistically significant differences among parameters. Metho dical parameter optimization will also be useful in this situation. 30 3. MODELLING OF URINARY BLADDER TISSUE 3.1 Introduction We hypothesize that mechanical properties of the urinary bladder differ between anatomical locations based on the findings of Korossis et al. (2009) and that the stress state of tissues samples is different in solutions of varying o smolarity. To test this theory, we performed uniaxial tests on porcine UB strips. Our data is in the form of 5 incremental stress relaxation curves. An outlin e of the experimental set up is first presented followed by each model that was tested. The bulk o f the analytical work for this data was in finding a suitable viscoelastic constitutive model to fit the data with. To make intentions for each model clear, t he methods, results, and discussion will be presented for each model before moving on to the next model. Four models were explored. The first two fit elastic and viscoelastic parameters: (1) Quasilinear viscoelasticity using a discrete relaxation spectrum (DRS) ( Figure 6 and Sections 1.3.3., 1.3.4. ) and (2) the Weichert Model ( Figure 4 and Section 1.3 .1.) using a continuous relaxation spectrum (CRS) ( Figure 6 and Section 1.3.4. ) in place of a series of infinite Maxwell elements. The last two models only fi t the viscoelastic parameters by fitting a reduced relaxation function (RRF) ( Eq. 32 and Section 1 .3.4. ) to normalized data. These models were (3) RRF using a CRS defined by Eq. 31 with n=1 and (4) RRF using a CRS defined by Eq. 31 with n=2. Advantages and disadvantages for each model will be discussed as well as future work that will improve on the re sults presented. 3.2. Experimental Methods The following experiment was performed by Tyler Tuttle, a fellow lab member of EMBR Lab and Mechanical Engineering PhD candidate at Michigan State University. Porcine urinary bladders were harvested for two different stress - relaxation experiments. Three bladders were used per each experiment (six total bladders harvested). Five samples of 1x3 cm strips cut from each bl adder, totaling in 15 samples per experiment. Samples were mounted on the uniaxial machine via a clamp and submerged in Krebs solution 18 hours prior to testing. 31 Figure 14. A UB sa mple mounted on the uniaxial machine, submerged in Krebs solution. The str ess - relaxation test consisted of an instantaneous strain of 0.25 followed by a wait period of 30 minutes, a second strain increase of 0.25 followed by a 45 minutes wait period, and three incremental strain increases of 0.5 held for 45 minutes each ( Figure 15 ). Figure 15. Five strains applied to the uniaxial machine (left) and the recorded stress relaxation behavior (right) The first experiment was a test to identify differen ces between anatomical locations of the viscoelastic properties of the UB . Samples were isolated from five locations: dorsal, ventral, lateral, lower body, and trigone ( Figure 16 ). The Krebs solution used was 300 mOsm/L for all tests . The s econd experiment was a test to identify the effect of varying Krebs solution osmolarity (i .e., swelling) on the viscoelastic behavior of the tissue . For this experiment, five samples were all isolated from the lateral location ( Figure 16 ). Four different o smolarity were tested, specifically 100, 300 (homeostatic) , 600, and 900 mOsm/L ; a fifth t est was performed dry (without Krebs solution). 32 Figure 16. Samples cut for the location mechanical test, specifically dorsal (D), ventral (V), lateral (L), lower b ody (LB), and trigone (T) locations (left) and s amples cut from the lateral location for t he swelling test (right) 3.3. Previous Work Prior to viscoelastic modelling of the data sets, analysis of the peak and relaxed stress was performed by Tyler Tuttle, who presented the following results at the 8 th World Congress of Biomechanics in 2018 [Tutt le and Roccabianca , 2018]. Average peak stresses and relaxed stresses are plotted for location trials ( Figure 17 ) and swelling trials ( Figure 18 ). A two - way ANOVA was performed to compare significant differences among the trials. The p - values for the resul ts are given for location ( Table 3 ) and for swelling ( Table 4 ). Figure 17. Average peak stresses for first three strains 0.25, 0.5, and 1.0 for peak stresses (left) and relaxed stresses (right) for location trials 33 Table 3. P - values comparing statisti cal significance between peak stresses (left) and relaxed stresses (right) based on a natomical location Figure 18. Average peak stresses for first three strains 0.25, 0.5, and 1.0 for peak stresses (left) and relaxed stresses (right) for swelling trials Table 4. P - values comparing statistical significance between peak stresses (left) and relaxed stresses (right) based o n bath osmolarity There were statistically significant differences (p<0.05) between dorsal/lower body, 900 and 300 mOsm/L, and dry and 300 mOsm/L for both peak and relaxed stresses. In general, the trigone and lower body behaved differently than the dorsal and ventral regions. Samples tested dry and with 900 mOsm/L had higher stresses than other swelling trials. These results motivate d us to fit parameters to the entire stress relaxation to see if viscoelastic behavior is different between anatomical region s and bath osmolarity. 34 3. 4 . Optimization Method Model parameters were fit to the data by minimizing the root mean square error (RMS E) between theoretical and experimental stress relaxation curves through an iterative process shown in Figure 1 9 . The minimization was performed using the MAT LA RMSE for Model (1) was based on the forces ( Eq. 36 ). Models (2) (4) c alc ulated RMSE based on stress ( Eq. 37 ). ( 36 ) ( 37 ) Figure 1 9 . Diagram of the optimization procedure 3.5. Model 1: Quasil inear Viscoelasticity in FEBi o, Single Stress - Relaxation 3.5.1. FEBio Model For this method, FEBio, an open source software developed by Musculoskeletal Research Laboratories at University of Utah, was used to create a finite element model of our experiment. The reasoning for using FE for a 1D problem (which does not require 3D modelling) was so that this work can translate to 3D models of the whole bladder in future projects. The material chosen was Quasilinear Viscoelasticity with a neo - Hookean e lastic component. The second Piola Kir chhoff stress for this material is defined as: 35 (38) where is the elastic stress component and G(t) is the relaxation functions defined as: (39) Note that this is not the same as a reduced relaxation function which follows the condition . The elastic copmonent is related to the strain energy function by Lastly, the strain energy function used here neo - Hookean , define d as : (40) w here the parameters are defined as ( Eq. 6 ), , , and . Our model used a relaxation function with n = 3, so the inputs to this model were , , . T he recorded experimental data was assumed to be the force acting on the thinnest cross - section of tissue ( Figure 20 ). The theoretical model forces at the cross section located on the Y - Z plane (surface B in Figure 21 and Figure 2 2 ) were calculated from the outputs of Y - displacement, Z - displacement, and X - stress . (41 ) For ce was used for the optimizat ion of the material parameters so that data could be compared directly from experiments , without conversion to stress based on changing cross - sectional area . 36 Figure 20 . Porcine UB mounted on the uniaxial machine (left) and an early FE model of the experiment showing stress in the x - direction (right) and the cross section where the forces were calculated as a dotted line Considering the symmetry of the problem, we included in t he FE model created in FEBio 1 /8 th of the original 1x1x0.5 cm sample ( Fig ure 21 ). A 7x7x7 elements mesh was used with bias in the x - direction to capture the stresses at the cross section of interest. The model assumes no slippage of the tissue in the clamps , as pictured above in Figure 20 , represented by a rigid contact between the tissue model and the rigid body (representing the clamp) . A displacement of the rigid body based on the first strain was then applied followed by 1800 seconds (30 minutes) of stress relaxation . Only the first strain was applied in FEBio. The parameter , and the relaxation strengths , associated with the relaxation times , from Eq. 39 . The relaxation times were constants, the reasoning for which was dis cussed in Section 2.5.3 . 37 Figure 21 . An isometric view of 1/8 th the original tissues with a 7x7x7 mesh (left) and the stress of the model immediately after the first strain where the color bar is stress [ Pa ] a nd surfaces A, B, and C are cross sections shown in Fig . 2 2 (right) Figure 2 2 . FEBio model during a relaxation phase showing 1/8 th of the original tissue at differe nt cross sections (labelled in Fig. 21 above) 38 3.5.2. Methods: Initial Results The f irst attempts to fit parameters to the FEBio model used the optimization procedure outlined in Figure 1 9 . U nfortunately, this method allowed the optimization to end prematu rely at local minima. We know this is a local minimum because manually fitting parameters was more precis e. An example of this is seen in Figure 23 . In this example, time steps assigned in FEBio were spaced 1 second apart for the first 100 seconds, and 10 seconds apart for the final 1700 seconds (for a total of 30 minutes). Figure 23. Optimization by procedure in Fig. 23 (RMSE = 6.29) versus hand fitted results (RMSE = 4.35) and raw data on a normal time s cale (left) and logarithmic time scale (right) The reason for MATLAB finding a false minimum became clear after performing a parametric study. Stress also decides the asymptote tha t stress relaxes to (elastic modulus times the strain) seen in Figure 15 . Stress is sensitive to any given relaxation pair only up until , one logarithmic decade following , as seen in Figure 2 4 . This is the time after when the relaxation pair are fully relaxed. The time marks the inflection point of the isolated r elaxation pair on a logarithmic time sca le. 39 Figure 24. Relaxation times (pictured , ) are the location of the inflection point, on the logarithmic time scale, of the individual exponentials (boxes) and they decay to approximately zero withi n one logarithmic decade (b ubbles ) eter at a time, increasing and decreasing its value, assessing the new function value (RMSE in this case), and altering the parameter once more based on the gradient it found for that single parameter. It repeats this process for every parameter until the function value has converged, meaning it has reached a local minimum. This does not work well for functions with multiple minima, as in the relaxation function we are using. The reason for this is that it does not consider locational parameter sensitivity (seen in Figure 13 ) , (2) inter - parameter sensitivity (discussed in Section 2.5.3. ), and (3) the magnitude at which it must change parameters to escape a local minimum. There may be many ways to approach this problem, such as a genetic algorithm (discussed in Section 2.4.3. ) that is equipped to handle parameter relationships. Another approach would be a weighted error function that gives equal importance to different time secti ons. For the initial optimization results, equally spaced than to say, which decays quickly after . The last approach would be to force the optimization process by fitting p ara meters only to the sections they have the largest impact on. We chose this approach over a weighted error function due to shorter run times of the code. 40 3.5.3. Methods: Guided Optimization for Quasil inear Viscoelasticity As discussed in the literature r evi ew, there may be infinite solutions to this model. Parameter optimization sectioning. To solve this issue, we implement ed a step - by - step optimizati on process that accounted parameter sensitivity to locations along the x - axis. RMSE is calculated for individual sections, most affected by certain parameters. To section the curve appropriately we used a geometric mean between time constants, which will b e denoted as . As seen in Figure 2 5 , the geometric mean between two times and is a good cutoff point as long as they are one logarithmic decade apart ( ). Figure 2 5 . The geometric mean between two time con stants is the logarithmic hallway point where the decay from is complete if . The steps of the optimization process were as follows: 1) 2) Fit to section between the and , where is the geometric mean defined as: 3) Repeat step (2) until reaching , fit from to 41 Figure 2 6 shows which step uses which section fo r a Prony Series with n = 3. must be chosen at least one logarithmic decade apart for the geometric mean to work as intended. Figure 26. Sections numbered in order that they are used for fitting parameters where section (1) is used for step (1), section (2) - (3) are used for step (2) where , and section (4) is used for step (3) These steps ensured that constants were fit only to portions of the curve that were sensitive to their changes . Time constants were chosen by manually f itting parameters. The time constants 2.5 s, 25 s, and 500 s were then used for every sample. 3.5.4. Results In the following table belong to 2.5 s, 25 s, and 500 s respectively, and E is - Hookean elastic component. 42 Table 5 . Average fitted parameters by anatomical location (top) and bath osmolarity (bottom) for the linear viscoelastic FEBio model (model 1) Table 6 . Average RMSE values for m odel (1) The dorsal location showed a high standard deviation for ameters had high standard deviations as well. The average parameters in Table 5 were used to plot the stress relaxation and RRF by rerunning the model through FEBio. The following figures show the average curves first for the location experiments (RRF in F igure 2 7 , stress - relaxation profile in Figure 2 9 ) and second for the swelling e xperiments (RRF in Figure 2 8 , stress - relaxation profile in Figure 30 ). 43 Figure 2 7 . Average RRF (G(T)) plotted for each anatomical location dorsal (D), ventral (V), lateral (L), lower body (B), and trigone (T) for model (1) Figure 2 8 . Average RRF plotted for each osmolarity trial 100, 300, 600, and 900 mOsm/L and dry for model (1) 44 Figure 2 9 . Average stress relaxation curves plotted for each anatomical location dorsal (D), ve ntral (V), lateral (L), lower body (B), and trigone (T) for model (1) 45 Figure 30. Average stress relaxation curves plotted for each osmolarity trial 100, 300, 600, and 900 mOsm/L and dry for model (1) 46 The RRF show different characteristics from the str ess relaxation curves. For example, the lateral location had the lowest pe ak and relaxed stresses but higher values in the RRF. A higher RRF is generally related with less stress - relaxation capability and decreased compliance. The dorsal location had the h ighest peak and relaxed stresses but had an average RRF compared to all lo cations. Samples with 100 and 900 mOsm/L had lower stress relaxation (higher RRF values), while at the same time these samples showed the lowest peak and relaxed stresses. 3.5.5. Dis cussio n One source of error for calculating the elastic component is that the tests might have been ramped to the next strain before reaching a fully relaxed state. In the experimental results, the most significant differences in peak and relaxed stresses by location and osmolarity were observed in the later pulls (i.e., at higher strains). This leads us to conclude that the parameters fitted to the first strain only might not reveal statistically significant differences between samples because of the low v alue of deformation applied. User bias is also introduced by using fixed relaxation times that might alter results with unknown magnitude of significance. Lastly, n oise from data collection is very high at the peak stresses and might skew results (as seen in higher standard deviation in the dorsal location). To reduce the number of parameters, reduce user bias, and fit data to all curves, a different approach was taken without FEBio. 3.6. Model 2: Weichert Model with Logarithmic Exponential Distributi on 3.6.1. Methods The strain history must be accounted for when modelling the incremental step strains of our stress - relaxation results. Therefore, the following model was derived from the Weichert model combined with a continuous relaxation function in pl ace of a series of Maxwell elements: ( 42 ) 47 ( 43 ) Whe re is the elastic modulus, is the elastic modulus of the viscoelastic component, is the current strain state, is the previous strain state, is stress as a function of time during the stress - relaxation st ep , is the stress as a function of time from the previous strain, and is a contin u ous distribution function. To minimize the number of parameters, a logarithmic exponential distribution function [ Tahmasbi and Rezaei , 2008 ] that uses tw o parameters ( was chosen to describe : ( 44 ) The parameters , , , and p were fit to the entire strain history via Eq. 43 and Eq. 44 . 3.6. 2 . Results Table 7 . Average fitted parameters fitted to Eq. 43 and Eq. 44 by anatomical location (top) and swelling (bottom ) for the logarithmic exponential CRS (model 2) 48 Table 8 . Average RMSE values for m odel ( 2 ) Figure 31 . The incremental stress - relaxation curves for location trials using parameters from Table 5 for model (2) Figure 3 2 . The incremental stress - relaxation curves for swelling trials using parameters from Table 5 for model (2) 49 The trends shown by the average parameters reflect what observed in the experimental results for osmolarity trials but not for location trials. The pe ak and relaxed stresses were highest for 900 and 100 mOsm/L and lowest for 300 and 600 mOsm/L which confirms what observed in the experim ental results. The lower body and trigone regions, however, had peak and relaxed stresses included within the range def ined by the other locations, while they showed to have the highest values in the experimental results. Additionally, the differences betw een parameters are not statistically significant. A two - way ANOVA was performed to analyze differences between all para meters, and none showed significance. This is no surprise, considering that some of the standard deviations for the parameter p are rathe r large (i.e., lower body, 100 mOsm/L, and 900 mOsm/L). Additionally, the RMSE values were rather high as can be seen i n Table 8 . 3.6.3. Discussion This model was found to be incapable of accurately fitting the data based on the high RMSE values and confli cting results for location trials . One limitation of this model is that it is not capable of capturing nonlinear elastic behavior of the elastic portion of the tissue . The theoretical curves increase in stress proportional to the strain (in the relaxed sta te) while the experimental stress seems to behave in a non - linear way . While some of the data were desc ribed rather accurately by this model, there were other datasets where the description was less accurate ( Figure 3 3 ). Figure 3 3 . Experimental versus t heoretical data for a good fit (bladder #2 900 mOsm/L) (left) and a bad fit (bladder #2 dry) (right) using model (2) 50 Other spectrums , such as the discrete spectrum for the Weichert model and gaussian CRS used by Nagatomi et al. (2008) , were combined with E q. 43 but had the same consequences as the ones seen in Figure 3 3 . Therefore, the Weichert model was abandoned for the remaining models while CRS continued to be tested . 3.7. Model 3: Reduced Relaxation Function with Log Normal Distribution ( ) 3.7.1. Methods There are two possible approaches to describe incremental strains in a stress - relaxation test on nonlinear viscoelastic materials employing an RRF, although this has been rarely done in literature (see Section 2.5.2. ). O ne approach co nsists of using QLV with but including a non - linear, hyperelastic component to describe the elastic portion of the response. The second approach is to model viscoelasticity by decoupling it from the nonlinear hyperelastic behavior. We will on ly use the sec ond approach for the remaining models and discuss how to implement the first approach in the section on future work. We normalize each stress relaxation curve (five per data set) and fit a CRS to all five simultaneously. Each stress relaxation curve was no rmalized by the following equation (from Sarver et al., 2003): (45) w here is the current stress - relaxation, is the final (relaxed) stress for the current curve, and is the peak stress for the current curve seen in Figure 3 4 . Figure 34 . Peak and relaxed stresses for each stress relaxation curve are used to normalize the curve by Eq. 45 [ Sarver et al., 2003 ] 51 The normalization technique allows us to study the re laxation behavior without modelling the hyperelastic component. Figure 3 5 shows how data normalized for all strains have similar relaxation behavior. Figure 3 5 . Normalized for all strains in the same experiment bladder #6 dorsal (left) and blad der #3, 100 mOsm/L (right) show no trends in increasing final stress that are seen in the nonnormalized data We employ here the dual gaussian relaxation spectrum with [ Nagatomi et al., 2008 ] and optimize the parameters : ( 31 ) ( 33 ) 52 3.7.2. Results Table 9 . Parameter s for Eq. 31 with n= 1 fit to data for location trials (top) and swelling trials (bottom) for m odel (3) Table 10 . A verage RMSE values for m odel (3) Figure 3 6 . Average relaxation behavior (RRF) by location based on parameters from Table 7 for m odel (3) 53 Figure 3 7 . Average relaxation behavior (RRF) by bath osmolarity based on parameters from Table 9 for m odel (3) 3.7. 3. Discussion T he logarithmic exponential distribution from M odel (2) was tested but had an asymptote at and could not be fit to the RRF. The CRS using a spectrum defined by Eq. 31 with matched predictions of behavior by location and osmola rity in terms of the overall RRF ( Fig. 3 6 and Fig. 3 7 ). Lower body and trigone regions had the lowest compliance among the location trials. The dry s amples had the lowest compliance, 100 and 900 mOsm/L samples were the next lowest compliant, and 300 and 60 0 mOsm/L samples had the highest compliance. These results seem to agree with what observed experimentally, however the parameters themselves are too large to be analyzed by statistical software. Amplitude is on the magnitude of . The values seem too large considering that the values reported in literature for were on the magnitude of [Nagatomi et al., 2008], so we suspect that this spec trum is not appropriate but that the normalization process improved results. 54 3. 8 . Model 4: Reduced Relaxation Function with Log Normal Distribution (n = 2) 3. 8 .1. Methods A second trial of th e dual gaussian spectrum ( Eq. 31 ) used and was able to fit results with smaller RMSE values on average than . Unfortunately, the double peaked spectrum is an ill - posed problem like the Prony Series, where there is no unique solution. A ny change to the mean , requires an adjustment of its dependent parameters and . A dditionally, the f itting process involves integration of the initial - gues s spectrum, making it an inverse problem very sensitive to noise in the data . An a is outlined in F igure 3 8 . While we are aware of the ill - p osedness of the problem , once the initial guess for each parameter is set, the solution is , if not unique, quite repeatable using this procedure. This assumption is supported by the empiric observation that reruns with th e same initial guess resulted in t he exact same results. For this reason, we worked on identifying . Initial values for the parameters were taken from Nagatomi et al. (2008) and are given below in Table 1 1 . Table 11 . Initial guess parameters from Nagatomi et al. (2008) Figure 3 8 . Guided optimization for Log - outlined in Fig. 1 9 55 3.8. 2. Results Table 12 . Parameters for Eq. 31 with n=2 fit to data for location trials (top) and bath osmolarity trials (bottom) for m odel (4) Table 13 . Average RMSE values for m odel (4) 56 Figure 3 9 . Average RRF based on parameters f rom Table 1 2 for location trials for m odel (4) Figure 40 . Average RRF based on parameters from Table 1 2 for swelling trials for m odel (4) 57 3.8.3. Discussion Results based on location strongly indicated lowered relaxation capability for lower body and trig one regions and 600 mOsm/L bath osmolarity. The results for bath osmolarity conflict with the results from the spectrum using . Although the location tria ls have plots that look like the expected behavior, we believe results are not representative because some of the standard deviations are rather large compared to the means. 3.9. Overall Conclusions There were observable tre nds in the relaxation behavior o f different locations. The trigone and lower body regions had lower stress relaxation (took longer to relax to the same stress) and higher final stresses. These trends are what we expect to see based on results from Korossis et al. for location. Similarly, swelling samples of the 900 mOsm/L, 100 mOsm/L, and dry tests had lowered stress relaxation and higher final stresses. While the trends are promising, the parameters that produce the average stress - relaxation curves failed p - tests. 3.9.1. Sources of Error One source of error for all trials is that strains were re - applied before the sample reached a fully relaxed state, only observable on the logarithmic scale. This can skew the resulting parameters of the elastic portion of each model. A source of error fo r method (1) would be that choosing the time constants through a manual fitting could influence the results for Figure 13 . To eliminate error from time constant selection, the interconversion algorithms mentioned in the literature review would provide the best time constants without user bia s because the discrete spectra are based on peaks in the continuous spectra ( Figure 12 ). Method (4) is insufficient to fit the parameters when the spectrum is bimodal ( ). The authors who used this model on the UB mention the use of a genetic algorithm to solve this [ Nagatomi et al. , 2004 ] . A genetic algorithm might be more effective in fit ting parameters for methods (1) and (4) because these the solutions are non - unique. 58 We speculate that an appropriate constitutive model might show statistically signif icant differences of viscoelastic parameters based on location and swelling. 3.9.2. Future Work and Works in Progress There are other time spectrums that were not tested for methods (3) and (4) that might better fit our data such as the straight value spec trum ( Eq. 46 ) which, when combined with equation Eq. 32 , gives a generalized reduced relaxed function in the form of Eq. 47 . ( 46 ) , . ( 47 ) A second approach to account for the nonlinear viscoelastic behavior would be to use the FEBio model from method (1) and replace the linear elastic component with a nonlinear component. Manually fitt ed trials using EFD Neo - Hookean Elasticity were run on th is model and show that it can fit relaxed stresses better than linear viscoelasticity. The downside is that the optimization process still arrives at false minima unsupervised (see Figure 41 ). This m odel will require the same guided optimization method use d for Model (1). Figure 41. Viscoelastic solid with EFD neo - Hookean elastic component optimized in MATLAB (left) and manually fitted (right) 59 The swelling trials can also be improved with multi - ph asic modelling. The FEBio model was re - designed as a trip hasic material with a solid mixture of Donnan Equilibrium and neo - Hookean elastic component. A fixed charge density of the solid matrix was ramped from zero to 100 [ ] during the first analys is step ( Fig. 42 . step A). The next step ramped the initial bath concentration of 300 mOsm/L from up to 900 mOsm/L for one run and down to 100 mOsm/L for a second run ( Fig. 42 . step B) . Both runs were then strained to 0.25 ( Fig. 42 . step C). The last step is where stress relaxation is observed from fluid moving in or out of the porous solid ( Fig. 42 . step C). This proof - of - concept model showed that the different bath concentration led to different peak stresses, like our experimental results ( Figure 43 ). Figure 42. FEBio triphasic model showing stress in the x - direction in different stages (A) loading the fixed charge density of the Donnan Equilibrium neo - Hookean solid (B) loading the bath osmolarity (C) applying a strain to the specimen (D ) the relaxation phase 60 Figure 43. The x - stress through the different stages of analysis of triphasic solid showing higher peak stress for ramping bath osmolarity from 300 to 900 mOsm/L (left) and lower peak stress for ramping bath osmolarity from 300 to 100 mOsm/L (right) The models presented in this paper showed a gradual improvement in the results and indicated the shortcoming of certain models. Combinations of different effective methods such as the normalization technique, DRS - CRS conversion, and t he use of nonl inear elastic models and triphasic models will aid in better understanding the relaxation behavior of the UB. Improved parameter characterization will lead to better full organ models in the future. 6 1 BIBLOGRAPHY 6 2 BIBLIOGRAPHY Abramowitch, Steven D., and Savio L - Y. Woo. "An improved method to analyze the stress relaxation of ligaments following a finite ramp time based on the quasi - linear viscoelastic theory." J. Biomech. Eng. 126.1 (2004): 92 - 97. Baumgaertel, M., and H. H. Winter. "Determination of discr ete rel axation and retardation time spectra from dynamic mechanical data." Rheologica Acta 28.6 (1989): 511 - 519. Chantereau, P., et al. "Mechanical properties of pelvic soft tissue of young women and impact of aging." International urogynecology journal 25 .11 (20 14): 1547 - 1553. Cost, Thomas L., and Eric B. Becker. "A multidata method of approximate Laplace transform inversion." International journal for numerical methods in engineering 2.2 (1970): 207 - 219. Dahms, S. E., et al. "Composition and biomechanica l prope rties of the bladder acellular matrix graft: comparative analysis in rat, pig and human." British journal of urology 82.3 (1998): 411 - 419. Doehring, Todd C., Evelyn O. Carew, and Ivan Vesely. "The effect of strain rate on the viscoelastic response o f aorti c valve tissue: a direct - fit approach." Annals of biomedical engineering 32.2 (2004): 223 - 232. Flack, Chandra, and C. R. Powell. "The worldwide economic impact of neurogenic bladder." Current bladder dysfunction reports 10.4 (2015): 350 - 354. Fry, A. Wagg, CH. "Visco - elastic properties of isolated detrusor smooth muscle." Scandinavian Journal of Urology and Nephrology 33.201 (1999): 12 - 18. Fung, Yuan - cheng. Biomechanics: mechanical properties of living tissues . Springer Science & Business Media, 1981 . Gupta, H. S., et al. "In situ multi - level analysis of viscoelastic deformation mechanisms in tendon collagen." Journal of structural biology 169.2 (2010): 183 - 191. Goh, S. M., M. N. Charalambides, and J. G. Williams. "Determination of the constitutive con stants of non - linear viscoelastic materials." Mechanics of Time - Dependent Materials 8.3 (2004): 255 - 268. Kohandel, M., S. Sivaloganathan, and Giuseppe Tenti. "Estimation of the quasi - linear viscoelastic parameters using a genetic algorithm." Mathematical a nd Comp uter Modelling 47.3 - 4 (2008): 266 - 270. Korossis, Sotirios, et al. "Regional biomechanical and histological characterisation of the passive porcine urinary bladder: Implications for augmentation and tissue engineering strategies." Biomaterials 30.2 ( 2009): 266 - 275. Labus, Kevin M., and Christian M. Puttlitz. "Viscoelasticity of brain corpus callosum in biaxial tension." Journal of the Mechanics and Physics of Solids 96 (2016): 591 - 604. Liu, Z., and K. Yeung. "The preconditioning and stress relaxation of skin tissue." Journal of Biomedical & Pharmaceutical Engineering 2.1 (2008): 22 - 28. Maas, Steve, et al. "FEBio theory manual." Musculoskeletal Research Laboratories, University of Utah, Salt Lake City, UT (2011). Martins, Pedro ALS, et al. "U niaxial mec hanical behavior of the human female bladder." International urogynecology journal 22.8 (2011): 991 - 995. 6 3 McDougall, Ian, Nese Orbey, and John M. Dealy. "Inferring meaningful relaxation spectra from experimental data." Journal of Rheology 58.3 (2 014): 779 - 7 97. Mun, Sungho, and Goangseup Zi. "Modeling the viscoelastic function of asphalt concrete using a spectrum method." Mechanics of Time - Dependent Materials 14.2 (2010): 191 - 202. Nagatomi, Jiro, et al. "Changes in the biaxial viscoelastic response of the uri nary bladder following spinal cord injury." Annals of biomedical engineering 32.10 (2004): 1409 - 1419. Nagatomi, Jiro, et al. "Contribution of the extracellular matrix to the viscoelastic behavior of the urinary bladder wall." Biomechanics and mo deling in m echanobiology 7.5 (2008): 395 - 404. Natali, A. N., et al. "Bladder tissue biomechanical behavior: experimental tests and constitutive formulation." Journal of biomechanics 48.12 (2015): 3088 - 3096. Park, S. W., and R. A. Schapery. "Methods of inte rconversion between linear viscoelastic material functions. Part I A numerical method based on Prony series." International journal of solids and structures 36.11 (1999): 1653 - 1675. Park, S. W., and Y. R. Kim. "Fitting Prony - series viscoelastic models with power - law presmoothing." Journal of materials in civil engineering 13.1 (2001): 26 - 32. Polliack, T., et al. "Clinical and economic consequences of volume - or time - dependent intermittent catheterization in patients with spinal cord lesions and neuropathic b ladder." Sp inal cord 43.10 (2005): 615 - 619. Roccabianca, Sara, and Tamara Reid Bush. "Understanding the mechanics of the bladder through experiments and theoretical models: Where we started and where we are heading." Technology 4.01 (2016): 30 - 41. Roylance , David. "E ngineering viscoelasticity." Department of Materials Science and Engineering Massachusetts Institute of Technology, Cambridge MA 2139 (2001): 1 - 37. Sarver, Joseph J., Paul S. Robinson, and Dawn M. Elliott. "Methods for quasi - linear viscoelastic modeling of soft tissue: application to incremental stress - relaxation experiments." J. Biomech. Eng. 125.5 (2003): 754 - 758. Schapery, Richard Allan. "A simple collocation method for fitting viscoelastic models to experimental data." (1962). Shanbhag, Sachi n. "pyReSpe ct: A computer program to extract discrete and continuous spectra from stress relaxation experiments ." Macromolecular Theory and Simulations 28.3 (2019): 1900005. Shen, Zhilei Liu, et al. "Viscoelastic properties of isolated collagen fibrils." Biophysical journal 100.12 (2011): 3008 - 3015. Tahmasbi, Rasool, and Sadegh Rezaei. "A two - parameter lifetime distribution with decreasing failure rate." Computational Statistics & Data Analysis 52.8 (2008): 3889 - 3901. Thiruchelvam , Nikesh, et al. "Neurotransmission an d viscoelasticity in the ovine fetal bladder after in utero bladder outflow obstruction." American Journal of Physiology - Regulatory, Integrative and Comparative Physiology 284.5 (2003): R1296 - R1305. Tschoegl, Nicholas W. The phenomenological theory of lin ear viscoelastic behavior: an introduction . Springer Science & Business Media, 1989 . 6 4 Tuttle, T . and Roccabianca , S Effects of swelling on urinary bladder wall mechanics th World Congress of Biomechanics, July 8 - 12, 2018, Convention Centre, Dublin, Ire land. Presentation. Van Mastrigt, Ron, and J. C. Nagtegaal. "Dependence of the viscoelastic response of the urinary bladder wall on strain rate." Medical and Biological Engineering and Computing 19.3 (1981): 291 - 296. Vincent, Julian. "Basic elasticity and viscoelasticity ." Structural biomaterials 1 (2012): 1 - 28. Zanetti, Elisabetta M., et al. "Bladder tissue passive response to monotonic and cyclic loading." Biorheology 49.1 (2012): 49 - 63.