Multiphasic finite element framework for modeling hydrated mixtures with multiple neutral and charged solutes

Modeling Fatigue Failure of Cartilage and Fibrous Biological Tissues Using Constrained Reactive Mixture Theory

Fatigue failure in biological soft tissues plays a critical role in the etiology of chronic soft tissue injuries and diseases such as osteoarthritis (OA). Understanding failure mechanisms is hindered by the decades-long timescales over which damage takes place. Analyzing the factors contributing to fatigue failure requires the help of validated computational models developed for soft tissues. This study presents a framework for fatigue failure of fibrous biological tissues based on reaction kinetics, where the composition of intact and fatigued material regions can evolve via degradation and breakage over time, in response to energy-based fatigue and damage criteria. Using reactive constrained mixture theory, material region mass fractions are governed by the axiom of mass balance. Progression of fatigue is controlled by an energy-based reaction rate, with user-selected probability functions defining the damage propensity of intact and fatigued material regions. Verification of this reactive theory, which is implemented in the open-source FEBio finite element software, is provided in this study. Validation is also demonstrated against experimental data, showing that predicted damage can be linked to results from biochemical assays. The framework is also applied to study fatigue failure during frictional contact of cartilage. Simulating previous experiments suggests that frictional effects slightly increase fatigue progression, but the main driver is cyclic compressive contact loading. This study demonstrated the ability of theoretical models to complement and extend experimental findings, advancing our understanding of the time progression of fatigue in biological tissues.

Simulating Cerebral Edema and Ischemia After Traumatic Acute Subdural Hematoma Using Triphasic Swelling Biomechanics

Poor outcome following traumatic acute subdural hematoma (ASDH) is associated with the severity of the primary injury and secondary injury including cerebral edema and ischemia. However, the underlying secondary injury mechanism contributing to elevated intracranial pressure (ICP) and high mortality rate remains unclear. Cerebral edema occurs in response to the exposure of the intracellular fixed charge density (FCD) after cell death, causing ICP to increase. The increased ICP from swollen tissue compresses blood vessels in adjacent tissue, restricting blood flow and leading to ischemic damage. We hypothesize that the mass occupying effect of ASDH exacerbates the ischemic injury, leading to ICP elevation, which is an indicator of high mortality rate in the clinic. Using FEBio (febio.org) and triphasic swelling biomechanics, this study modeled clinically relevant ASDHs and simulated post-traumatic brain swelling and ischemia to predict ICP. Results showed that common convexity ASDH significantly increased ICP by exacerbating ischemic injury, and surgical removal of the convexity ASDH may control ICP by preventing ischemia progression. However, in cases where the primary injury is very severe, surgical intervention alone may not effectively decrease ICP, as the contribution of the hematoma to the elevated ICP is insignificant. In addition, interhemispheric ASDH, located between the cerebral hemispheres, does not significantly exacerbate ischemia, supporting the conservative surgical management generally recommended for interhemispheric ASDH. The joint effect of the mass occupying effect of the blood clot and resulting ischemia contributes to elevated ICP which may increase mortality. Our novel approach may improve the fidelity of predicting patient outcome after motor vehicle crashes and traumatic brain injuries due to other causes.

Region-Dependent Mechanical Properties of Human Brain Tissue Under Large Deformations Using Inverse Finite Element Modeling

This study aims to facilitate intracranial simulation of traumatic events by determining the mechanical properties of different anatomical structures of the brain. Our experimental indentation paradigm used fresh, post-operative human tissue, which is highly advantageous in determining mechanical properties without being affected by postmortem time. This study employed an inverse finite element approach coupled with experimental indentation data to characterize mechanical properties of the human hippocampus (CA1, CA3, dentate gyrus), cortex white matter, and cortex grey matter. We determined that an uncoupled viscoelastic Ogden constitutive formulation was most appropriate to represent the mechanical behavior of these different regions of brain. Anatomical regions were significantly different in their mechanical properties. The cortex white matter was stiffer than cortex grey matter, and the CA1 and dentate gyrus were both stiffer than cortex grey matter. Although no sex dependency was observed, there were trends indicating that male brain regions were generally stiffer than corresponding female regions. In addition, there were no statistically significant age dependent differences. This study provides a structure-specific description of fresh human brain tissue mechanical properties, which will be an important step toward explicitly modeling the heterogeneity of brain tissue deformation during TBI through finite element modeling.

Finite Element Implementation of Computational Fluid Dynamics With Reactive Neutral and Charged Solute Transport in FEBio

The objective of this study was to implement a novel fluid-solutes solver into the open-source finite element software FEBio, that extended available modeling capabilities for biological fluids and fluid-solute mixtures. Using a reactive mixture framework, this solver accommodates diffusion, convection, chemical reactions, electrical charge effects, and external body forces, without requiring stabilization methods that were deemed necessary in previous computational implementations of the convection-diffusion-reaction equation at high Peclet numbers. Verification and validation problems demonstrated the ability of this solver to produce solutions for Peclet numbers as high as 1011, spanning the range of physiological conditions for convection-dominated solute transport. This outcome was facilitated by the use of a formulation that accommodates realistic values for solvent compressibility, and by expressing the solute mass balance such that it properly captured convective transport by the solvent and produced a natural boundary condition of zero diffusive solute flux at outflow boundaries. Since this numerical scheme was not necessarily foolproof, guidelines were included to achieve better outcomes that minimize or eliminate the potential occurrence of numerical artifacts. The fluid-solutes solver presented in this study represents an important and novel advancement in the modeling capabilities for biomechanics and biophysics as it allows modeling of mechanobiological processes via the incorporation of chemical reactions involving neutral or charged solutes within dynamic fluid flow. The incorporation of charged solutes in a reactive framework represents a significant novelty of this solver. This framework also applies to a broader range of nonbiological applications.

Modeling Inelastic Responses Using Constrained Reactive Mixtures

This study reviews the progression of our research, from modeling growth theories for cartilage tissue engineering, to the formulation of constrained reactive mixture theories to model inelastic responses in any solid material, such as theories for damage mechanics, viscoelasticity, plasticity, and elasto-plastic damage. In this framework, multiple solid generations α can co-exist at any given time in the mixture. The oldest generation is denoted by α = s and is called the master generation, whose reference configuration X s is observable. The solid generations α are all constrained to share the same velocity v s , but may have distinct reference configurations X α . An important element of this formulation is that the time-invariant mapping F α s = ∂ X α / ∂ X s between these reference configurations is a function of state, whose mathematical formulation is postulated by constitutive assumption. Thus, reference configurations X α are not observable ( α ≠ s ) . This formulation employs only observable state variables, such as the deformation gradient F s of the master generation and the referential mass concentrations ρ r α of each generation, in contrast to classical formulations of inelastic responses which rely on internal state variable theory, requiring evolution equations for those hidden variables. In constrained reactive mixtures, the evolution of the mass concentrations is governed by the axiom of mass balance, using constitutive models for the mass supply densities ρ ˆ r α . Classical and constrained reactive mixture approaches share considerable mathematical analogies, as they both introduce a multiplicative decomposition of the deformation gradient, also requiring evolution equations to track some of the state variables. However, they also differ at a fundamental level, since one adopts only observable state variables while the other introduces hidden state variables. In summary, this review presents an alternative foundational approach to the modeling of inelastic responses in solids, grounded in the classical framework of mixture theory.

Cartilage biomechanics: From the basic facts to the challenges of tissue engineering

Articular cartilage (AC) is the thin tissue that covers the long bone ends in the joints and that ensures the transmission of forces between adjacent bones while allowing nearly frictionless movements between them. AC repair is a technologic and scientific challenge that has been addressed with numerous approaches. A major deadlock is the capacity to take in account its complex mechanical properties in repair strategies. In this review, we first describe the major mechanical behaviors of AC for the non-specialists. Then, we show how researchers have progressively identified specific mechanical parameters using mathematical models. There are still gaps in our understanding of some of the observations concerning AC biomechanical properties, particularly the differences in extracellular matrix stiffness measured at the microscale and at the millimetric scale. Nevertheless, for bioengineering applications, AC repair strategies must take into account what are commonly considered the main mechanical features of cartilage: its ability to withstand high stresses through three main behaviors (elasticity, poroelasticity and swelling). Finally, we emphasize that future studies need to investigate AC mechanical properties at different scales, particularly the gradient of mechanical properties around cells and across the cartilage depth, and the differences in mechanical properties at different scales. This multi-scale approach could greatly enhance the success of AC restorative approaches.

Generation and evaluation of input values for computational analysis of transport processes within tissue cultures

Techniques for tissue culture have seen significant advances during the last decades and novel 3D cell culture systems have become available. To control their high complexity, experimental techniques and their Digital Twins (modelling and computational tools) are combined to link different variables to process conditions and critical process parameters. This allows a rapid evaluation of the expected product quality. However, the use of mathematical simulation and Digital Twins is critically dependent on the precise description of the problem and correct input parameters. Errors here can lead to dramatically wrong conclusions. The intention of this review is to provide an overview of the state-of-the-art and remaining challenges with respect to generating input values for computational analysis of mass and momentum transport processes within tissue cultures. It gives an overview on relevant aspects of transport processes in tissue cultures as well as modelling and computational tools to tackle these problems. Further focus is on techniques used for the determination of cell-specific parameters and characterization of culture systems, including sensors for on-line determination of relevant parameters. In conclusion, tissue culture techniques are well-established, and modelling tools are technically mature. New sensor technologies are on the way, especially for organ chips. The greatest remaining challenge seems to be the proper addressing and handling of input parameters required for mathematical models. Following Good Modelling Practice approaches when setting up and validating computational models is, therefore, essential to get to better estimations of the interesting complex processes inside organotypic tissue cultures in the future.

The intrinsic quality of proteoglycans, but not collagen fibres, degrades in osteoarthritic cartilage

The function of articular cartilage as a load-bearing connective tissue is derived primarily from a balanced interaction between the swelling proteoglycan (PG) matrix and tension-resistant collagen fibrous network. Such balance is compromised during joint disease such as osteoarthritis (OA) due to degradation to PGs and/or collagens. While the PG degradation is generally thought to be related to a loss of protein abundance, the collagenous degradation is more complex as it can be caused independently by a decrease of collagen content, disorganisation of fibrous structure and softening of individual collagen fibrils. A comprehensive understanding of the initial trajectories of degradation of PGs and collagen network can improve our chance of finding potential therapeutic solutions for OA. Here, we developed geometrically, structurally, and compositionally realistic and sample-specific Finite Element (FE) models under the framework of multiphasic mixture theory, from which the elastic moduli of collagen fibres and the PG load-bearing quality in healthy and diseased cartilages were estimated by numerical optimisation of the multi-step indentation stress relaxation force-time curves. We found the intrinsic quality of collagen fibres, measured by their elastic moduli, to stay constant for healthy and diseased cartilages. Combining with previous findings which show unaltered collagen content during early stages of OA, our results suggest the disorganisation of collagen fibrous network as the first form of collagenous degradation in osteoarthritic cartilage. We also found that PG degradation involves not only a loss of protein abundance, but also the quality of the remaining PGs in generating sufficient osmotic pressure for load bearing. This study sheds light on the mechanism of OA pathogenesis and highlights the restoration of collageneous organisation in cartilage as key medical intervention for OA. STATEMENT OF SIGNIFICANCE: Collagen network in articular cartilage consists of individual fibres that are organised into depth-dependent structure specialised for joint load-bearing and lubrication. During osteoarthritis, the collagen network undergoes mechanical degradation, but it is unclear if a loss of content, disorganisation of fibrous structure, or softening of individual fibres causes this degeneration. Using mechanical indentation, Finite Element modelling, and numerical optimisation methods, we determined that individual fibres did not soften in early disease stage. Together with previous findings showing unaltered collagen content, our results pinpoint the disorganisation of collagen structure as the main culprit for early collagenous degradation in osteoarthritic cartilage. Thus, early restoration in cartilage of collagen organisation, instead of individual fibre quality, may be key to slow osteoarthritis development.

Modelling articular cartilage: the relative motion of two adjacent poroviscoelastic layers

In skeletal joints two layers of adjacent cartilage are often in relative motion. The individual cartilage layers are often modelled as a poroviscoelastic material. To model the relative motion, noting the separation of scales between the pore level and the macroscale, a homogenization based on multiple scale asymptotic analysis has been used in this study to derive a macroscale model for the relative translation of two poroviscoelastic layers separated by a very thin layer of fluid. In particular the fluid layer thickness is essentially zero at the macroscale so that the two poroviscoelastic layers are effectively in contact and their interaction is captured in the derived model via a set of interfacial conditions, including a generalization of the Beavers-Joseph condition at the interface between a viscous fluid and a porous medium. In the simplifying context of a uniform geometry, constant fixed charge density, a Newtonian interstitial fluid and a viscoelastic scaffold, modelled via finite deformation theory, we present preliminary simulations that may be used to highlight predictions for how oscillatory relative movement of cartilage under load influences the peak force the cartilage experiences and the extent of the associated deformations. In addition to highlighting such cartilage mechanics, the systematic derivation of the macroscale models will enable the study of how nanoscale cartilage physics, such as the swelling pressure induced by fixed charges, manifests in cartilage mechanics at much higher lengthscales.

A Hybrid Biphasic Mixture Formulation for Modeling Dynamics in Porous Deformable Biological Tissues

The primary aim of this study is to establish the theoretical foundations for a solid-fluid biphasic mixture domain that can accommodate inertial effects and a viscous interstitial fluid, which can interface with a dynamic viscous fluid domain. Most mixture formulations consist of constituents that are either all intrinsically incompressible or compressible, thereby introducing inherent limitations. In particular, mixtures with intrinsically incompressible constituents can only model wave propagation in the porous solid matrix, whereas those with compressible constituents require internal variables, and related evolution equations, to distinguish the compressibility of the solid and fluid under hydrostatic pressure. In this study, we propose a hybrid framework for a biphasic mixture where the skeleton of the porous solid is intrinsically incompressible but the interstitial fluid is compressible. We define a state variable as a measure of the fluid volumetric strain. Within an isothermal framework, the Clausius-Duhem inequality shows that a function of state arises for the fluid pressure as a function of this strain measure. We derive jump conditions across hybrid biphasic interfaces, which are suitable for modeling hydrated biological tissues. We then illustrate this framework using confined compression and dilatational wave propagation analyses. The governing equations for this hybrid biphasic framework reduce to those of the classical biphasic theory whenever the bulk modulus of the fluid is set to infinity and inertia terms and viscous fluid effects are neglected. The availability of this novel framework facilitates the implementation of finite element solvers for fluid-structure interactions at interfaces between viscous fluids and porous-deformable biphasic domains, which can include fluid exchanges across those interfaces.

On the use of constrained reactive mixtures of solids to model finite deformation isothermal elastoplasticity and elastoplastic damage mechanics

This study presents a framework for plasticity and elastoplastic damage mechanics by treating materials as reactive solids whose internal composition evolves in response to applied loading. Using the framework of constrained reactive mixtures, plastic deformation is accounted for by allowing loaded bonds within the material to break and reform in a stressed state. Bonds which break and reform represent a new generation with a new reference configuration, which is time-invariant and provided by constitutive assumption. The constitutive relation for the reference configuration of each generation may depend on the selection of a suitable yield measure. The choice of this measure and the resulting plastic flow conditions are constrained by the Clausius-Duhem inequality. We show that this framework remains consistent with classical plasticity approaches and principles. Verification of this reactive plasticity framework, which is implemented in the open source FEBio finite element software (febio.org), is performed against standard 2D and 3D benchmark problems. Damage is incorporated into this reactive framework by allowing loaded bonds to break permanently according to a suitable damage measure, where broken bonds can no longer store free energy. Validation is also demonstrated against experimental data for problems involving plasticity and plastic damage. This study demonstrates that it is possible to formulate simple elastoplasticity and elastoplastic damage models within a consistent framework which uses measures of material mass composition as theoretically observable state variables. This theoretical frame can be expanded in scope to account for more complex behaviors.

Finite Element Implementation of Biphasic-Fluid Structure Interactions in febio

In biomechanics, solid-fluid mixtures have commonly been used to model the response of hydrated biological tissues. In cartilage mechanics, this type of mixture, where the fluid and solid constituents are both assumed to be intrinsically incompressible, is often called a biphasic material. Various physiological processes involve the interaction of a viscous fluid with a porous-hydrated tissue, as encountered in synovial joint lubrication, cardiovascular mechanics, and respiratory mechanics. The objective of this study was to implement a finite element solver in the open-source software febio that models dynamic interactions between a viscous fluid and a biphasic domain, accommodating finite deformations of both domains as well as fluid exchanges between them. For compatibility with our recent implementation of solvers for computational fluid dynamics (CFD) and fluid-structure interactions (FSI), where the fluid is slightly compressible, this study employs a novel hybrid biphasic formulation where the porous skeleton is intrinsically incompressible but the fluid is also slightly compressible. The resulting biphasic-FSI (BFSI) implementation is verified against published analytical and numerical benchmark problems, as well as novel analytical solutions derived for the purposes of this study. An illustration of this BFSI solver is presented for two-dimensional (2D) airflow through a simulated face mask under five cycles of breathing, showing that masks significantly reduce air dispersion compared to the no-mask control analysis. In addition, we model three-dimensional (3D) blood flow in a bifurcated carotid artery assuming porous arterial walls and verify that mass is conserved across all fluid-permeable boundaries. The successful formulation and implementation of this BFSI solver offers enhanced multiphysics modeling capabilities that are accessible via an open-source software platform.

Direct Osmotic Pressure Measurements in Articular Cartilage Demonstrate Nonideal and Concentration-Dependent Phenomena

The osmotic pressure in articular cartilage serves an important mechanical function in healthy tissue. Its magnitude is thought to play a role in advancing osteoarthritis. The aims of this study were to: (1) isolate and quantify the magnitude of cartilage swelling pressure in situ; and (2) identify the effect of salt concentration on material parameters. Confined compression stress-relaxation testing was performed on 18 immature bovine and six mature human cartilage samples in solutions of varying osmolarities. Direct measurements of osmotic pressure revealed nonideal and concentration-dependent osmotic behavior, with magnitudes approximately 1/3 those predicted by ideal Donnan law. A modified Donnan constitutive behavior was able to capture the aggregate behavior of all samples with a single adjustable parameter. Results of curve-fitting transient stress-relaxation data with triphasic theory in febio demonstrated concentration-dependent material properties. The aggregate modulus HA increased threefold as the external concentration decreased from hypertonic 2 M to hypotonic 0.001 M NaCl (bovine: HA=0.420±0.109 MPa to 1.266±0.438 MPa; human: HA=0.499±0.208 MPa to 1.597±0.455 MPa), within a triphasic theory inclusive of osmotic effects. This study provides a novel and simple analytical model for cartilage osmotic pressure which may be used in computational simulations, validated with direct in situ measurements. A key finding is the simultaneous existence of Donnan osmotic and Poisson-Boltzmann electrostatic interactions within cartilage.

Size-Dependent Effective Diffusivity in Healthy Human and Porcine Joint Synovium

Intra-articular drug delivery can be effective in targeting a diseased joint but is hampered by rapid clearance times from the diarthrodial joint. The synovium is a multi-layered tissue that surrounds the diarthrodial joint and governs molecular transport into and out of the joint. No models of drug clearance through synovium exist to quantify diffusivity across solutes, tissue type and disease pathology. We previously have developed a finite element model of synovium as a porous, permeable, fluid-filled tissue and used an inverse method to determine urea's effective diffusivity (D_eff) in de-vitalized synovium explants.²² Here we apply this method to determine D_eff from unsteady diffusive transport of model solutes and confirm the role of molecular weight in solute transport. As molecular weight increased, D_eff decreased in both human and porcine tissues, with similar behavior across the two species. Unsteady transport was well-described by a single exponential transient decay in concentration, yielding solute half-lives (t_1/2) that compared favorably with the D_eff determined from the finite element model fit. Determined values for D_eff parallel prior observations of size-dependent in vivo drug clearance and provide an intrinsic parameter with greater ability to resolve size-dependence in vitro. Thus, this work forms the basis for understanding the influence of size on drug transport in synovium and can guide future studies to elucidate the role of charge and tissue pathology on the transport of therapeutics in healthy and pathological human synovium.

Finite deformation elastography of articular cartilage and biomaterials based on imaging and topology optimization

Tissues and engineered biomaterials exhibit exquisite local variation in stiffness that defines their function. Conventional elastography quantifies stiffness in soft (e.g. brain, liver) tissue, but robust quantification in stiff (e.g. musculoskeletal) tissues is challenging due to dissipation of high frequency shear waves. We describe new development of finite deformation elastography that utilizes magnetic resonance imaging of low frequency, physiological-level (large magnitude) displacements, coupled to an iterative topology optimization routine to investigate stiffness heterogeneity, including spatial gradients and inclusions. We reconstruct 2D and 3D stiffness distributions in bilayer agarose hydrogels and silicon materials that exhibit heterogeneous displacement/strain responses. We map stiffness in porcine and sheep articular cartilage deep within the bony articular joint space in situ for the first time. Elevated cartilage stiffness localized to the superficial zone is further related to collagen fiber compaction and loss of water content during cyclic loading, as assessed by independent T2 measurements. We additionally describe technical challenges needed to achieve in vivo elastography measurements. Our results introduce new functional imaging biomarkers, which can be assessed nondestructively, with clinical potential to diagnose and track progression of disease in early stages, including osteoarthritis or tissue degeneration.

Combined Experimental Approach and Finite Element Modeling of Small Molecule Transport Through Joint Synovium to Measure Effective Diffusivity

Trans-synovial solute transport plays a critical role in the clearance of intra-articularly (IA) delivered drugs. In this study, we present a computational finite element model (FEM) of solute transport through the synovium validated by experiments on synovial explants. Unsteady diffusion of urea, a small uncharged molecule, was measured through devitalized porcine and human synovium using custom-built diffusion chambers. A multiphasic computational model was constructed and optimized with the experimental data to extract effective diffusivity for urea within the synovium. A monotonic decrease in urea concentration was observed in the donor bath over time, with an effective diffusivity found to be an order of magnitude lower in synovium versus that measured in free solution. Parametric studies incorporating an intimal cell layer with varying thickness and varying effective diffusivities were performed, revealing a dependence of drug clearance kinetics on both parameters. The findings of this study indicate that the synovial matrix impedes urea solute transport out of the joint with little retention of the solute in the matrix.

Simulating cerebral edema and delayed fatality after traumatic brain injury using triphasic swelling biomechanics

Objectives: Contemporary finite element (FE) models, like that from the Global Human Body Models Consortium (GHBMC), have been useful for developing safety systems to reduce the severity of injuries in motor vehicle crashes (MVCs), including traumatic brain injury (TBI). However, not all injury occurs during the MVC. Cerebral edema after TBI contributes to mortality by increasing intracranial pressure (ICP) and preventing adequate cerebral blood supply. The focus of this study was to model post-traumatic cerebral edema and subsequent mortality due to increased ICP.Methods: Brain tissue swells in a manner consistent with triphasic biomechanics, which models biological tissues as a charged deformable porous solid matrix (fixed charge density [FCD]), a solvent, and monovalent counter-ions (cerebrospinal fluid). Fluid uptake into the brain is driven by the Gibbs-Donnan osmotic pressure as the FCD is exposed when cells die. Post-TBI edema was simulated in FEBio (febio.org), which includes triphasic material formulations.The GHBMC mesh was imported into FEBio, and each element was assigned a FCD to represent impact-related cell death based on its maximum principal strain (MPS) experienced during the crash-simulation using the stock GHBMC model and LS-DYNA. The ensuing pathophysiology was simulated in FEBio in two steps. First, the brain swelled in response to exposure of FCD, causing some adjacent elements to compress as fluid was redistributed. Biologically, the compression was assumed to reduce blood flow and cause ischemic cell death, represented by additional exposure of FCD, swelling, and increased ICP. Using published prognostic models of clinical outcome, mortality was predicted based on ICP.Results: Post-traumatic volume ratio of elements ranged from less than 30% (compaction) to greater than 200% (swelling). Predicted ICP values for a fatal impact were as high as 8.55 kPa (64.1 mmHg), which is associated with a 99% probability of death.Conclusion: To the best of our knowledge, this is the first study to simulate post-traumatic brain swelling to predict outcome. By incorporating swelling, ischemia, and cell death, our novel approach may improve fidelity of predicting outcome after MVCs. A strength of our approach is relying on the validated GHBMC model to predict brain deformation in the crash-scenario. The main goal of the current study was to demonstrate feasibility of simulating post-injury swelling using triphasic biomechanics. We successfully predicted clinically relevant increases in ICP that suggest a high likelihood of death when simulating a fatal impact scenario, however, more validation of our methodology is needed.

Finite Element Framework for Computational Fluid Dynamics in FEBio

The mechanics of biological fluids is an important topic in biomechanics, often requiring the use of computational tools to analyze problems with realistic geometries and material properties. This study describes the formulation and implementation of a finite element framework for computational fluid dynamics (CFD) in FEBio, a free software designed to meet the computational needs of the biomechanics and biophysics communities. This formulation models nearly incompressible flow with a compressible isothermal formulation that uses a physically realistic value for the fluid bulk modulus. It employs fluid velocity and dilatation as essential variables: The virtual work integral enforces the balance of linear momentum and the kinematic constraint between fluid velocity and dilatation, while fluid density varies with dilatation as prescribed by the axiom of mass balance. Using this approach, equal-order interpolations may be used for both essential variables over each element, contrary to traditional mixed formulations that must explicitly satisfy the inf-sup condition. The formulation accommodates Newtonian and non-Newtonian viscous responses as well as inviscid fluids. The efficiency of numerical solutions is enhanced using Broyden's quasi-Newton method. The results of finite element simulations were verified using well-documented benchmark problems as well as comparisons with other free and commercial codes. These analyses demonstrated that the novel formulation introduced in FEBio could successfully reproduce the results of other codes. The analogy between this CFD formulation and standard finite element formulations for solid mechanics makes it suitable for future extension to fluid-structure interactions (FSIs).

Analyzing the effects of instillation volume on intravesical delivery using biphasic solute transport in a deformable geometry

Ailments of the bladder are often treated via intravesical delivery-direct application of therapeutic into the bladder through a catheter. This technique is employed hundreds of thousands of times every year, but protocol development has largely been limited to empirical determination. Furthermore, the numerical analyses of intravesical delivery performed to date have been restricted to static geometries and have not accounted for bladder deformation. This study uses a finite element analysis approach with biphasic solute transport to investigate several parameters pertinent to intravesical delivery including solute concentration, solute transport properties and instillation volume. The volume of instillation was found to have a substantial impact on the exposure of solute to the deeper muscle layers of the bladder, which are typically more difficult to reach. Indeed, increasing the instillation volume from 50-100 ml raised the muscle solute exposure as a percentage of overall bladder exposure from 60-70% with higher levels achieved for larger instillation volumes. Similar increases were not seen for changes in solute concentration or solute transport properties. These results indicate the role that instillation volume may play in targeting particular layers of the bladder during an intravesical delivery.

A Formulation for Fluid Structure-Interactions in FEBio Using Mixture Theory

Many physiological systems involve strong interactions between fluids and solids, posing a signicant challenge when modeling biomechanics. The objective of this study was to implement a fluid-structure interaction (FSI) solver in the free, open-source finite element code FEBio (febio.org), that combined the existing solid mechanics and rigid body dynamics solver with a recently-developed computational fluid dynamics (CFD) solver. A novel Galerkin-based finite element FSI formulation was introduced based on mixture theory, where the FSI domain was described as a mixture of fluid and solid constituents that have distinct motions. The mesh was defined on the solid domain, specialized to have zero mass, negligible stiffness and zero frictional interactions with the fluid, whereas the fluid was modeled as isothermal and compressible. The mixture framework provided the foundation for evaluating material time derivatives in a material frame for the solid and in a spatial frame for the fluid. Similar to our recently reported CFD solver, our FSI formulation did not require stabilization methods to achieve good convergence, producing a compact set of equations and code implementation. The code was successfully verified against benchmark problems and an analytical solution for squeeze-film lubrication. It was validated against experimental measurements of the flow rate in a peristaltic pump, and illustrated using non-Newtonian blood flow through a bifurcated carotid artery with a thick arterial wall. The successful formulation and implementation of this FSI solver enhances the multiphysics modeling capabilities in FEBio relevant to the biomechanics and biophysics communities.

Intracellular ion and protein nanoparticle-induced osmotic pressure modify astrocyte swelling and brain edema in response to glutamate stimuli

Intracellular tension activity plays a crucial role in cytotoxic brain edema and astrocyte swelling. Here, a few genetically encoded FRET-based tension probes were designed to detect cytoskeletal structural tension optically, including their magnitude and vectors. The astrocyte swelling resulted in GFAP tension increment, which is associated with the antagonistic effect of inward microfilaments (MFs) and microtubules (MTs) forces. In glutamate-induced astrocyte swelling, GFAP tension rise resulted from outward ion and protein nanoparticle-induced osmotic pressure (PN-OP) increases, where PN-OP could be elicited by MF and MT depolymerization, protein nanoparticle production, and activation of cofilin and stathmin-1. Attenuation of both ion osmotic pressure and PN-OP by drug combinations, together with free-radical scavenger, relieved cerebral edema in vivo. The study suggests that intracellular osmotic pressure (especially PN-OP) has a pivotal role in glutamate-induced astrocyte swelling and brain edema. Recovery of cytoplasmic potential is a promising target to develop new drugs and cure brain edema.

Prediction of probability of fatality due to brain injury in traffic accidents

Objective: Fatal brain injuries result from physiological changes in brain tissues, subsequent to primary damage caused by head impact. Although efforts have been made in past studies to estimate the probability of brain injury, none of them involved prediction of such physiological changes. The goal of this study was to evaluate the fatality prediction capability of a novel approach that predicts an increase in intracranial pressure (ICP) due to primary head injury to estimate the fatality rate using clinical data that correlate ICP with fatality rate. Methods: A total of 12 sets of head acceleration time histories were used to represent no, severe, and fatal brain injury. They were obtained from the literature presenting head kinematics data in noninjurious volunteer sled tests or from accident reconstruction for severe and fatal injury cases. These were first applied to a Global Human Body Models Consortium (GHBMC) head-brain model to predict nodal displacement time histories of the brain, which were then fed into FEBio to predict ICP. A Weibull distribution was applied to the data for the relationship between fatality rate and ICP obtained from a clinical paper to estimate fatality rate from ICP (procedure A). Fatality rate was also estimated by applying the temporal and spatial maximum value of maximum principal strain (MPS_max) obtained from the GHBMC simulation to an injury probability function for MPS_max (procedure B). Estimated fatality rates were compared between the 2 procedures. Results: Both procedures estimated higher average fatality rate for higher injury severity. The average fatality rate for procedure A without ischemia representation and procedure B was 72.4 and 51.0% for the fatal injury group and 8.2 and 21.7% for the severe injury group, respectively, showing that procedure A provides more distinct classification between fatal and nonfatal brain injury. It was also found that representation of ischemia in procedure A provides results sensitive to injury severity and impact conditions, requiring further validation of the initial estimate for the relationship between brain compression and ischemic cell death. Conclusions: Prediction of the probability of fatality by means of a combination of simulations of the primary brain deformation and subsequent ICP increase was found to be more distinct compared to the prediction of primary injury alone combined with the injury probability function from a past study in the select 12 head impact cases.

How can 50 years of solute transport data in articular cartilage inform the design of arthritis therapeutics?

OBJECTIVE

For the last half century, transport of nutrients and therapeutics in articular cartilage has been studied with various in vitro systems that attempt to model in vivo conditions. However, experimental technique, tissue species, and tissue storage condition (fresh/frozen) vary widely and there is debate on the most appropriate model system. Additionally, there is still no clear overarching framework with which to predict solute transport properties based on molecular characteristics. This review aims to develop such a framework, and to assess whether experimental procedure affects trends in transport data.

METHODS

Solute data from 31 published papers that investigated transport in healthy articular cartilage were obtained and analyzed for trends.

RESULTS

Here, we show that diffusivity of spherical and globular solutes in cartilage can be predicted by molecular weight (MW) and hydrodynamic radius via a power-law relationship. This relationship is robust for many solutes, spanning 5 orders of magnitude in MW and was not affected by variations in cartilage species, age, condition (fresh/frozen), and experimental technique. Traditional models of transport in porous media exhibited mixed effectiveness at predicting diffusivity in cartilage, but were good in predicting solute partition coefficient.

CONCLUSION

Ultimately, these robust relationships can be used to accurately predict and improve transport of solutes in adult human cartilage and enable the development of better optimized arthritis therapeutics.

Numerical Study on Electromechanics in Cartilage Tissue with Respect to Its Electrical Properties

Hyaline cartilage undergoes many substantial age-related physiochemical and biomechanical changes that reduce its ability to overcome the effects of mechanical stress and injury. In quest of therapeutic options, magnetic stimulation and electrical stimulation (ES) have been proposed for improving tissue engineering approaches for the repair of articular cartilage. The aim of this study is to summarize in silico investigations involving induced electrical properties of cartilage tissue due to various biophysical stimuli along their respective mathematical descriptions. Based on these, a preliminary numerical study involving electromechanical transduction in bovine cartilage tissue has been carried out using an open source finite element computational software. The simulation results have been compared to experimental results from the literature. This study serves as a basis for further in silico studies to better understand the behavior of hyaline cartilage tissue due to ES and to find an optimal stimulation protocol for the cartilage regeneration. Moreover, it provides an overview of the basic models along with mathematical description and scope for future research regarding electrical behavior of the cartilage tissue using open source software.

A Plugin Framework for Extending the Simulation Capabilities of FEBio

The FEBio software suite is a set of software tools for nonlinear finite element analysis in biomechanics and biophysics. FEBio employs mixture theory to account for the multiconstituent nature of biological materials, integrating the field equations for irreversible thermodynamics, solid mechanics, fluid mechanics, mass transport with reactive species, and electrokinetics. This communication describes the development and application of a new "plugin" framework for FEBio. Plugins are dynamically linked libraries that allow users to add new features and to couple FEBio with other domain-specific software applications without modifying the source code directly. The governing equations and simulation capabilities of FEBio are reviewed. The implementation, structure, use, and application of the plugin framework are detailed. Several example plugins are described in detail to illustrate how plugins enrich, extend, and leverage existing capabilities in FEBio, including applications to deformable image registration, constitutive modeling of biological tissues, coupling to an external software package that simulates angiogenesis using a discrete computational model, and a nonlinear reaction-diffusion solver. The plugin feature facilitates dissemination of new simulation methods, reproduction of published results, and coupling of FEBio with other domain-specific simulation approaches such as compartmental modeling, agent-based modeling, and rigid-body dynamics. We anticipate that the new plugin framework will greatly expand the range of applications for the FEBio software suite and thus its impact.

Cytoskeletal Contribution to Cell Stiffness Due to Osmotic Swelling; Extending the Donnan Equilibrium

Cell volume regulation is commonly analyzed with a model of a closed semipermeable membrane filled with impermeant mobile solutes and the Donnan Equilibrium is used to predict the hydrostatic pressure. This traditional model ignores the fact that most cells are filled with a crosslinked cytoskeleton that is elastic and can be stretched or compressed like a sponge with no obvious need to move mobile solutes. However, calculations show that under osmotic stress, the elastic energy of the cytoskeleton is far greater than the elastic energy of the membrane. Here we expand the traditional Donnan model to include the elasticity of a cytoskeleton with fixed charges and show that cell stiffening happens without a membrane.

Finite Element Formulation of Multiphasic Shell Elements for Cell Mechanics Analyses in FEBio

With the recent implementation of multiphasic materials in the open-source finite element (FE) software FEBio (febio.org), 3D models of cells embedded within the tissue may now be analyzed, accounting for porous solid matrix deformation, transport of interstitial fluid and solutes, membrane potential, and reactions. The cell membrane is a critical component in cell models, which selectively regulates the transport of fluid and solutes in the presence of large concentration and electric potential gradients, while also facilitating the transport of various proteins. The cell membrane is much thinner than the cell; therefore, in an FE environment, shell elements formulated as 2D surfaces in 3D space would be preferred for modeling the cell membrane, for the convenience of mesh generation from image-based data, especially for convoluted membranes. However, multiphasic shell elements are yet to be developed in the FE literature and commercial FE software. This study presents a novel formulation of multiphasic shell elements and its implementation in FEBio. The shell model includes front- and back-face nodal degrees of freedom for the solid displacement, effective fluid pressure and effective solute concentrations, and a linear interpolation of these variables across the shell thickness. This formulation was verified against classical models of cell physiology and validated against reported experimental measurements in chondrocytes. This implementation of passive transport of fluid and solutes across multiphasic membranes makes it possible to model the biomechanics of isolated cells or cells embedded in their extracellular matrix, accounting for solvent and solute transport.

Simulation of surface strain in tibiofemoral cartilage during walking for the prediction of collagen fiber orientation

The collagen fibers in the superficial layer of tibiofemoral articular cartilage exhibit distinct patterns in orientation revealed by split lines. In this study, we introduce a simulation framework to predict cartilage surface loading during walking to investigate if split line orientations correspond with principal strain directions in the cartilage surface. The two-step framework uses a multibody musculoskeletal model to predict tibiofemoral kinematics which are then imposed on a deformable surface model to predict surface strains. The deformable surface model uses absolute nodal coordinate formulation (ANCF) shell elements to represent the articular surface and a system of spring-dampers and internal pressure to represent the underlying cartilage. Simulations were performed to predict surface strains due to osmotic pressure, loading induced by walking, and the combination of both loading due to pressure and walking. Time-averaged magnitude-weighted first principal strain directions agreed well with split line maps from the literature for both the osmotic pressure and combined cases. This result suggests there is indeed a connection between collagen fiber orientation and mechanical loading, and indicates the importance of accounting for the pre-strain in the cartilage surface due to osmotic pressure.

Multiscale mechanics of the cervical facet capsular ligament, with particular emphasis on anomalous fiber realignment prior to tissue failure

The facet capsular ligaments encapsulate the bilateral spinal facet joints and are common sources of painful injury due to afferent innervation. These ligaments exhibit architectural complexity, which is suspected to contribute to the experimentally observed lack of co-localization between macroscopic strain and microstructural tissue damage. The heterogeneous and multiscale nature of this ligament, combined with challenges in experimentally measuring its microscale mechanics, hinders the ability to understand sensory mechanisms under normal or injurious loading. Therefore, image-based, subject-specific, multiscale finite-element models were constructed to predict the mechanical responses of the human cervical facet capsular ligament under uniaxial tensile stretch. The models precisely simulated the force-displacement responses for all samples ([Formula: see text]) and showed promise in predicting the magnitude and location of peak regional strains at two different displacements. Yet, there was a loss of agreement between the model and experiment in terms of fiber organization at large tissue stretch, possibly due to a lack of accounting for tissue failure. The mean fiber stretch ratio predicted by the models was found to be significantly higher in regions that exhibited anomalous fiber realignment experimentally than in regions with normal realignment ([Formula: see text]). The development of microstructural abnormalities was associated with the predicted fiber-level stretch ([Formula: see text]), but not with the elemental maximum principal stress or maximum principal strain by logistic regression. The multiscale models elucidate a potential mechanical basis for predicting injury-prone tissue domains and for defining the relationships between macroscopic ligament stretch and microscale pathophysiology in the subfailure regime.

FEBio: History and Advances

The principal goal of the FEBio project is to provide an advanced finite element tool for the biomechanics and biophysics communities that allows researchers to model mechanics, transport, and electrokinetic phenomena for biological systems accurately and efficiently. In addition, because FEBio is geared toward the research community, the code is designed such that new features can be added easily, thus making it an ideal tool for testing novel computational methods. Finally, because the success of a code is determined by its user base, integral goals of the FEBio project have been to offer support and outreach to our community; to provide mechanisms for dissemination of results, models, and data; and to encourage interaction between users. This review presents the history of the FEBio project, from its initial developments through its current funding period. We also present a glimpse into the future of FEBio.

Cyclic Mechanical Loading Enhances Transport of Antibodies Into Articular Cartilage

The goal of this study was to characterize antibody penetration through cartilage tissue under mechanical loading. Mechanical stimulation aids in the penetration of some proteins, but this effect has not characterized molecules such as antibodies (>100 kDa), which may hold some clinical value for treating osteoarthritis (OA). For each experiment, fresh articular cartilage plugs were obtained and exposed to fluorescently labeled antibodies while under cyclic mechanical load in unconfined compression for several hours. Penetration of these antibodies was quantified using confocal microscopy, and finite element (FE) simulations were conducted to predict fluid flow patterns within loaded samples. Transport enhancement followed a linear trend with strain amplitude (0.25-5%) and a nonlinear trend with frequency (0.25-2.60 Hz), with maximum enhancement found to be at 5% cyclic strain and 1 Hz, respectively. Regions of highest enhancement of transport within the tissue were associated with the regions of highest interstitial fluid velocity, as predicted from finite-element simulations. Overall, cyclic compression-enhanced antibody transport by twofold to threefold. To our knowledge, this is the first study to test how mechanical stimulation affects the diffusion of antibodies in cartilage and suggest further study into other important factors regarding macromolecular transport.

An Experimental and Finite Element Protocol to Investigate the Transport of Neutral and Charged Solutes across Articular Cartilage

Osteoarthritis (OA) is a debilitating disease that is associated with degeneration of articular cartilage and subchondral bone. Degeneration of articular cartilage impairs its load-bearing function substantially as it experiences tremendous chemical degradation, i.e. proteoglycan loss and collagen fibril disruption. One promising way to investigate chemical damage mechanisms during OA is to expose the cartilage specimens to an external solute and monitor the diffusion of the molecules. The degree of cartilage damage (i.e. concentration and configuration of essential macromolecules) is associated with collisional energy loss of external solutes while moving across articular cartilage creates different diffusion characteristics compared to healthy cartilage. In this study, we introduce a protocol, which consists of several steps and is based on previously developed experimental micro-Computed Tomography (micro-CT) and finite element modeling. The transport of charged and uncharged iodinated molecules is first recorded using micro-CT, which is followed by applying biphasic-solute and multiphasic finite element models to obtain diffusion coefficients and fixed charge densities across cartilage zones.

Swelling of Collagen-Hyaluronic Acid Co-Gels: An In Vitro Residual Stress Model

Tissue-equivalents (TEs), simple model tissues with tunable properties, have been used to explore many features of biological soft tissues. Absent in most formulations however, is the residual stress that arises due to interactions among components with different unloaded levels of stress, which has an important functional role in many biological tissues. To create a pre-stressed model system, co-gels were fabricated from a combination of hyaluronic acid (HA) and reconstituted Type-I collagen (Col). When placed in solutions of varying osmolarity, HA-Col co-gels swell as the HA imbibes water, which in turn stretches (and stresses) the collagen network. In this way, co-gels with residual stress (i.e., collagen fibers in tension and HA in compression) were fabricated. When the three gel types tested here were immersed in hypotonic solutions, pure HA gels swelled the most, followed by HA-Col co-gels; no swelling was observed in pure collagen gels. The greatest swelling rates and swelling ratios occurred in the lowest salt concentration solutions. Tension on the collagen component of HA-Col co-gels was calculated from a stress balance and increased nonlinearly as swelling increased. The swelling experiment results were in good agreement with the stress predicted by a fibril network + non-fibrillar interstitial matrix computational model.

Kinetics of charged antibiotic penetration into human intervertebral discs: A numerical study

Little quantitative information exists on the kinetics of charged antibiotic penetration into human intervertebral discs (IVD). This information is crucial for determining the dosage to use, timing of administration, and duration of treatment for infected IVDs. The objective of this study was to quantitatively analyze the transport of various charged antibiotics into human lumbar IVDs. Penetration of charged and uncharged antibiotics into a human lumbar disc was analyzed using a 3D finite element model. The valence (z) of the electrical charge of antibiotics varied from z=+2 (positively charged) to z=-2 (negatively charged). An uncharged antibiotic (z=0) was used as a control. Cases with intravenous (IV) administrations of different charged antibiotics were simulated. Our results showed that the electrical charge had great effects on kinetics of an antibiotic penetration into the IVD; with higher concentrations and uptakes for positively charged antibiotics than those for negatively charged ones. This study provides quantitative information on selecting antibiotics for treating intervertebral disc infections.

An overview of multiphase cartilage mechanical modelling and its role in understanding function and pathology

There is a long history of mathematical and computational modelling with the objective of understanding the mechanisms governing cartilage׳s remarkable mechanical performance. Nonetheless, despite sophisticated modelling development, simulations of cartilage have consistently lagged behind structural knowledge and thus the relationship between structure and function in cartilage is not fully understood. However, in the most recent generation of studies, there is an emerging confluence between our structural knowledge and the structure represented in cartilage modelling. This raises the prospect of further refinement in our understanding of cartilage function and also the initiation of an engineering-level understanding for how structural degradation and ageing relates to cartilage dysfunction and pathology, as well as informing the potential design of prospective interventions. Aimed at researchers entering the field of cartilage modelling, we thus review the basic principles of cartilage models, discussing the underlying physics and assumptions in relatively simple settings, whilst presenting the derivation of relatively parsimonious multiphase cartilage models consistent with our discussions. We proceed to consider modern developments that start aligning the structure captured in the models with observed complexities. This emphasises the challenges associated with constitutive relations, boundary conditions, parameter estimation and validation in cartilage modelling programmes. Consequently, we further detail how both experimental interrogations and modelling developments can be utilised to investigate and reduce such difficulties before summarising how cartilage modelling initiatives may improve our understanding of cartilage ageing, pathology and intervention.

Parameter-dependent behavior of articular cartilage: 3D mechano-electrochemical computational model

BACKGROUND AND OBJECTIVE

Changes in mechano-electrochemical properties of articular cartilage play an essential role in the majority of cartilage diseases. Despite of this importance, the specific effect of each parameter into tissue behavior remains still obscure. Parametric computational modeling of cartilage can provide some insights into this matter, specifically the study of mechano-electrochemical properties variation and their correlation with tissue swelling, water and ion fluxes. Thus, the aim of this study is to evaluate the influence of the main mechanical and electrochemical parameters on the determination of articular cartilage behavior by a parametric analysis through a 3D finite element model.

METHODS

For this purpose, a previous 3D mechano-electrochemical model, developed by the same authors, of articular cartilage behavior has been used. Young's modulus, Poisson coefficient, ion diffusivities and ion activity coefficients variations have been analyzed and quantified through monitoring tissue simulated response.

RESULTS

Simulation results show how Young's modulus and Poisson coefficient control tissue behavior rather than electrochemical properties. Meanwhile, ion diffusivity and ion activity coefficients appear to be vital in controlling velocity of incoming and outgoing fluxes.

CONCLUSIONS

This parametric study establishes a basic guide when defining the main properties that are essential to be included into computational modeling of articular cartilage providing a helpful tool in tissue simulations.

Multiscale Mechanical Model of the Pacinian Corpuscle Shows Depth and Anisotropy Contribute to the Receptor's Characteristic Response to Indentation

Cutaneous mechanoreceptors transduce different tactile stimuli into neural signals that produce distinct sensations of touch. The Pacinian corpuscle (PC), a cutaneous mechanoreceptor located deep within the dermis of the skin, detects high frequency vibrations that occur within its large receptive field. The PC is comprised of lamellae that surround the nerve fiber at its core. We hypothesized that a layered, anisotropic structure, embedded deep within the skin, would produce the nonlinear strain transmission and low spatial sensitivity characteristic of the PC. A multiscale finite-element model was used to model the equilibrium response of the PC to indentation. The first simulation considered an isolated PC with fiber networks aligned with the PC's surface. The PC was subjected to a 10 μm indentation by a 250 μm diameter indenter. The multiscale model captured the nonlinear strain transmission through the PC, predicting decreased compressive strain with proximity to the receptor's core, as seen experimentally by others. The second set of simulations considered a single PC embedded epidermally (shallow) or dermally (deep) to model the PC's location within the skin. The embedded models were subjected to 10 μm indentations at a series of locations on the surface of the skin. Strain along the long axis of the PC was calculated after indentation to simulate stretch along the nerve fiber at the center of the PC. Receptive fields for the epidermis and dermis models were constructed by mapping the long-axis strain after indentation at each point on the surface of the skin mesh. The dermis model resulted in a larger receptive field, as the calculated strain showed less indenter location dependence than in the epidermis model.

In Silico Investigation of Angiogenesis with Growth and Stress Generation Coupled to Local Extracellular Matrix Density

Mechanical interactions during angiogenesis, i.e., traction applied by neovessels to the extracellular matrix and the corresponding deformation, are important regulators of growth and neovascularization. We have previously designed, implemented, and validated a coupled model of angiogenesis in which a discrete microvessel growth model interacts with a continuous finite element mesh through the application of local remodeling sprout stresses (Edgar et al. in Biomech Model Mechanobiol, 2014). However, the initial implementation of this framework does not take matrix density into account when determined these remodeling stresses and is therefore insufficient for the study of angiogenesis within heterogeneous matrix environments such as those found in vivo. The objective of this study was to implement sensitivity to matrix density in the active stress generation within AngioFE in order to allow the study of angiogenic growth within a heterogeneous density environment. We accomplished this by scaling active sprout stresses relative to local matrix density using a scaling factor previously determined from experimental data. We then exercised the new functionality of the model by simulating angiogenesis within four different scenarios: homogeneous density, a narrow gap model, and matrix density gradient, and a construct subjected to repeated loading/unloading and preconditioning. These numerical experiments predicted heterogeneous matrix density in the initially homogeneous case, the closure and alignment of microvessels along a low-density gap, the formation of a unique cap-like structure during angiogenesis within a density gradient, and the alignment of microvessels in the absence of applied load due to preconditioning. The result of these in silico investigations demonstrate how matrix heterogeneity affects neovascularization and matrix deformation and provides a platform for studying angiogenesis in complicated and multi-faceted mechanical environments that microvessels experience in vivo.

Multiscale cartilage biomechanics: technical challenges in realizing a high-throughput modelling and simulation workflow

Understanding the mechanical environment of articular cartilage and chondrocytes is of the utmost importance in evaluating tissue damage which is often related to failure of the fibre architecture and mechanical injury to the cells. This knowledge also has significant implications for understanding the mechanobiological response in healthy and diseased cartilage and can drive the development of intervention strategies, ranging from the design of tissue-engineered constructs to the establishment of rehabilitation protocols. Spanning multiple spatial scales, a wide range of biomechanical factors dictate this mechanical environment. Computational modelling and simulation provide descriptive and predictive tools to identify multiscale interactions, and can lead towards a greater comprehension of healthy and diseased cartilage function, possibly in an individualized manner. Cartilage and chondrocyte mechanics can be examined in silico, through post-processing or feed-forward approaches. First, joint-tissue level simulations, typically using the finite-element method, solve boundary value problems representing the joint articulation and underlying tissue, which can differentiate the role of compartmental joint loading in cartilage contact mechanics and macroscale cartilage field mechanics. Subsequently, tissue-cell scale simulations, driven by the macroscale cartilage mechanical field information, can predict chondrocyte deformation metrics along with the mechanics of the surrounding pericellular and extracellular matrices. A high-throughput modelling and simulation framework is necessary to develop models representative of regional and population-wide variations in cartilage and chondrocyte anatomy and mechanical properties, and to conduct large-scale analysis accommodating a multitude of loading scenarios. However, realization of such a framework is a daunting task, with technical difficulties hindering the processes of model development, scale coupling, simulation and interpretation of the results. This study aims to summarize various strategies to address the technical challenges of post-processing-based simulations of cartilage and chondrocyte mechanics with the ultimate goal of establishing the foundations of a high-throughput multiscale analysis framework. At the joint-tissue scale, rapid development of regional models of articular contact is possible by automating the process of generating parametric representations of cartilage boundaries and depth-dependent zonal delineation with associated constitutive relationships. At the tissue-cell scale, models descriptive of multicellular and fibrillar architecture of cartilage zones can also be generated in an automated fashion. Through post-processing, scripts can extract biphasic mechanical metrics at a desired point in the cartilage to assign loading and boundary conditions to models at the lower spatial scale of cells. Cell deformation metrics can be extracted from simulation results to provide a simplified description of individual chondrocyte responses. Simulations at the tissue-cell scale can be parallelized owing to the loosely coupled nature of the feed-forward approach. Verification studies illustrated the necessity of a second-order data passing scheme between scales and evaluated the role that the microscale representative volume size plays in appropriately predicting the mechanical response of the chondrocytes. The tools summarized in this study collectively provide a framework for high-throughput exploration of cartilage biomechanics, which includes minimally supervised model generation, and prediction of multiscale biomechanical metrics across a range of spatial scales, from joint regions and cartilage zones, down to that of the chondrocytes.

Toward patient-specific articular contact mechanics

The mechanics of contacting cartilage layers is fundamentally important to understanding the development, homeostasis and pathology of diarthrodial joints. Because of the highly nonlinear nature of both the materials and the contact problem itself, numerical methods such as the finite element method are typically incorporated to obtain solutions. Over the course of five decades, we have moved from an initial qualitative understanding of articular cartilage material behavior to the ability to perform complex, three-dimensional contact analysis, including multiphasic material representations. This history includes the development of analytical and computational contact analysis methods that now provide the ability to perform highly nonlinear analyses. Numerical implementations of contact analysis based on the finite element method are rapidly advancing and will soon enable patient-specific analysis of joint contact mechanics using models based on medical image data. In addition to contact stress on the articular surfaces, these techniques can predict variations in strain and strain through the cartilage layers, providing the basis to predict damage and failure. This opens up exciting areas for future research and application to patient-specific diagnosis and treatment planning applied to a variety of pathologies that affect joint function and cartilage homeostasis.

A coupled model of neovessel growth and matrix mechanics describes and predicts angiogenesis in vitro

During angiogenesis, sprouting microvessels interact with the extracellular matrix (ECM) by degrading and reorganizing the matrix, applying traction forces, and producing deformation. Morphometric features of the resulting microvascular network are affected by the interaction between the matrix and angiogenic microvessels. The objective of this study was to develop a continuous-discrete modeling approach to simulate mechanical interactions between growing neovessels and the deformation of the matrix in vitro. This was accomplished by coupling an existing angiogenesis growth model which uses properties of the ECM to regulate angiogenic growth with the nonlinear finite element software FEBio (www.febio.org). FEBio solves for the deformation and remodeling of the matrix caused by active stress generated by neovessel sprouts, and this deformation was used to update the ECM into the current configuration. After mesh resolution and parameter sensitivity studies, the model was used to accurately predict vascular alignment for various matrix boundary conditions. Alignment primarily arises passively as microvessels convect with the deformation of the matrix, but active alignment along collagen fibrils plays a role as well. Predictions of alignment were most sensitive to the range over which active stresses were applied and the viscoelastic time constant in the material model. The computational framework provides a flexible platform for interpreting in vitro investigations of vessel-matrix interactions, predicting new experiments, and simulating conditions that are outside current experimental capabilities.

Computational modeling of chemical reactions and interstitial growth and remodeling involving charged solutes and solid-bound molecules

Mechanobiological processes are rooted in mechanics and chemistry, and such processes may be modeled in a framework that couples their governing equations starting from fundamental principles. In many biological applications, the reactants and products of chemical reactions may be electrically charged, and these charge effects may produce driving forces and constraints that significantly influence outcomes. In this study, a novel formulation and computational implementation are presented for modeling chemical reactions in biological tissues that involve charged solutes and solid-bound molecules within a deformable porous hydrated solid matrix, coupling mechanics with chemistry while accounting for electric charges. The deposition or removal of solid-bound molecules contributes to the growth and remodeling of the solid matrix; in particular, volumetric growth may be driven by Donnan osmotic swelling, resulting from charged molecular species fixed to the solid matrix. This formulation incorporates the state of strain as a state variable in the production rate of chemical reactions, explicitly tying chemistry with mechanics for the purpose of modeling mechanobiology. To achieve these objectives, this treatment identifies the specific theoretical and computational challenges faced in modeling complex systems of interacting neutral and charged constituents while accommodating any number of simultaneous reactions where reactants and products may be modeled explicitly or implicitly. Several finite element verification problems are shown to agree with closed-form analytical solutions. An illustrative tissue engineering analysis demonstrates tissue growth and swelling resulting from the deposition of chondroitin sulfate, a charged solid-bound molecular species. This implementation is released in the open-source program FEBio ( www.febio.org ). The availability of this framework may be particularly beneficial to optimizing tissue engineering culture systems by examining the influence of nutrient availability on the evolution of inhomogeneous tissue composition and mechanical properties, the evolution of construct dimensions with growth, the influence of solute and solid matrix electric charge on the transport of cytokines, the influence of binding kinetics on transport, the influence of loading on binding kinetics, and the differential growth response to dynamically loaded versus free-swelling culture conditions.

Synthesis rates and binding kinetics of matrix products in engineered cartilage constructs using chondrocyte-seeded agarose gels

Large-sized cartilage constructs suffer from inhomogeneous extracellular matrix deposition due to insufficient nutrient availability. Computational models of nutrient consumption and tissue growth can be utilized as an efficient alternative to experimental trials to optimize the culture of large constructs; models require system-specific growth and consumption parameters. To inform models of the [bovine chondrocyte]-[agarose gel] system, total synthesis rate (matrix accumulation rate+matrix release rate) and matrix retention fractions of glycosaminoglycans (GAG), collagen, and cartilage oligomeric matrix protein (COMP) were measured either in the presence (continuous or transient) or absence of TGF-β3 supplementation. TGF-β3's influences on pyridinoline content and mechanical properties were also measured. Reversible binding kinetic parameters were characterized using computational models. Based on our recent nutrient supplementation work, we measured glucose consumption and critical glucose concentration for tissue growth to computationally simulate the culture of a human patella-sized tissue construct, reproducing the experiment of Hung et al. (2003). Transient TGF-β3 produced the highest GAG synthesis rate, highest GAG retention ratio, and the highest binding affinity; collagen synthesis was elevated in TGF-β3 supplementation groups over control, with the highest binding affinity observed in the transient supplementation group; both COMP synthesis and retention were lower than those for GAG and collagen. These results informed the modeling of GAG deposition within a large patella construct; this computational example was similar to the previous experimental results without further adjustments to modeling parameters. These results suggest that these nutrient consumption and matrix synthesis models are an attractive alternative for optimizing the culture of large-sized constructs.