Multiscale Kinematic Growth Coupled With Mechanosensitive Systems Biology in Open-Source Software.

Multiscale Kinematic Growth Coupled With Mechanosensitive Systems Biology in Open-Source Software - PMC Skip to main content An official website of the United States government Here's how you know Here's how you know Search Log in Primary site navigation Search Log in Search PMC Full-Text Archive Search in PMC Journal List User Guide As a library, NLM provides access to scientific literature. Inclusion in an NLM database does not imply endorsement of, or agreement with, the contents by NLM or the National Institutes of Health. Learn more: PMC Disclaimer | PMC Copyright Notice J Biomech Eng . 2025 Apr 28;147(6):061004. doi: 10.1115/1.4068290 Multiscale Kinematic Growth Coupled With Mechanosensitive Systems Biology in Open-Source Software Steven A LaBelle 1, 6, #, ✉ , Mohammadreza Soltany Sadrabadi 2, 7, #, ✉ , Seungik Baek 3, 8 , Mohammad R K Mofrad 4 , Jeffrey A Weiss 1, 2 , Amirhossein Arzani 5, 3 Author information Article notes Copyright and License information 1 Department of Biomedical Engineering, University of Utah, Salt Lake City, UT 84112; Scientific Computing and Imaging Institute, University of Utah, Salt Lake City, UT 84112 2 Department of Mechanical Engineering, Northern Arizona University, Flagstaff, AZ 86011 3 Department of Mechanical Engineering, Michigan State University, East Lansing, MI 48824 4 Department of Bioengineering and Mechanical Engineering, University of California, Berkeley, Berkeley, CA 94720 5 Scientific Computing and Imaging Institute, University of Utah, Salt Lake City, UT 84112; Department of Mechanical Engineering, University of Utah, Salt Lake City, UT 84112 6 University of Utah 7 Northern Arizona University 8 Michigan State University SB3C conference. 1 S. A. LaBelle and M. Soltany Sadrabadi contributed equally to this work. 2 Corresponding author. e-mail: jeff.weiss@utah.edu 3 Corresponding author. e-mail: amir.arzani@sci.utah.edu ✉ Corresponding author. Received 2024 Sep 10; Revised 2025 Mar 14; Issue date 2025 Jun 1. Copyright © 2025 by ASME 0148-0731/2025/147(6)/061004/19/$25.00 PMC Copyright notice PMCID: PMC12147932  PMID: 40131342 Abstract Multiscale coupling between cell-scale biology and tissue-scale mechanics is a promising approach for modeling disease growth. In such models, tissue-level growth and remodeling (G&R) are driven by cell-level signaling pathways and systems biology models, where each model operates at different scales. Herein, we generate multiscale G&R models to capture the associated multiscale connections. At the cell-scale, we consider systems biology models in the form of systems of ordinary differential equations (ODEs) and partial differential equations (PDEs) representing the reactions between the biochemicals causing the growth based on mass-action or logic-based Hill-type kinetics. At the tissue-scale, we employ kinematic growth in continuum frameworks. Two illustrative test problems (a tissue graft and aneurysm growth) are examined with various chemical signaling networks, boundary conditions, and mechano-chemical coupling strategies. We extend two open-source software frameworks— febio and fenics —to disseminate examples of multiscale growth and remodeling simulations. One-way and two-way coupling between the systems biology and the growth models are compared and the effect of biochemical diffusivity and ODE versus PDE-based systems biology modeling on the G&R results are studied. The results show that growth patterns emerge from reactions between biochemicals, the choice between ODEs and PDEs systems biology modeling, and the coupling strategy. Cross-verification confirms that results for febio and fenics are nearly identical. We hope that these open-source tools will support reproducibility and education within the biomechanics community. Keywords: cell signaling, systems mechanobiology, growth and remodeling, multiscale modeling, aneurysm 1 Introduction Growth and remodeling (G&R) is the continuous process whereby biological tissues undergo changes in mass, geometry, mechanical properties, and function [ 1 , 2 ]. This inherently multiscale and multiphysics process is fundamental to biological phenomena such as morphogenesis, wound healing, homeostasis, and disease progression. G&R may be induced by biological processes (e.g., inflammation), mechanical signaling (e.g., tension, compression, shear, fixation), and systemic adaptations (e.g., in response to age, pregnancy, disease, smoking, or (epi-)genetics [ 3 , 4 ]). Biochemical and mechanical stimuli span physical scales ranging from subcellular systems to organs, and temporal scales ranging from milliseconds to decades [ 5 ]. Computational approaches have been developed to predict the G&R of various tissues in response to spatiotemporal variations in biochemical and mechanical stimuli with varying levels of biochemomechanical coupling. Recently, systems biology approaches have been used to couple cell-level dynamics with tissue-level adaptations. Such approaches have the potential to improve the understanding of diseases and identify therapeutic candidates based on patient specific factors such as demographics, genetics, or comorbidities. Despite these developments, numerous barriers preclude the integration of systems approaches to multidisciplinary research groups and educators in biomedical disciplines. Modeling is generally performed using proprietary software or in-house code, with sparse sharing of code or model input files. Modeling environments either lack user-friendly interfaces or rely on the user to compile code or adapt scripts in niche coding languages such as Fortran (e.g., abaqus subroutines). This is further complicated by the multiscale, multiphysical, and problem-specific nature of G&R, which necessitates tedious linearization and discretization of nonlinear partial differential equations (PDEs). Thus, our objective is to provide user-friendly, open-source software for simulating G&R. We cross-verify our implementation among two open-source software packages ( febio and fenics ). We demonstrate common concepts in G&R with simplified example problems that highlight the effects of various coupling formulations between biological and mechanical signals. Our efforts are outlined as follows: The remainder of the introduction provides a brief review of the relevant literature surrounding computational G&R as well as a review of relevant open-source software. The methods first formalize the kinematics, mechanics, and reactive transport. Next, software implementations, test problems, and cross-verification strategies are defined. The results overview findings from the test problems. The discussion considers the implications and limitations of the test problems and software implementations. Finally, the appendix provides useful mathematical relationships and detailed methods on the nonlinear finite element methods employed. Differential growth–the variation of growth rate and orientation between neighboring tissue subvolumes–emerges from cell-level phenomena, creating spatially heterogeneous tissue-level growth [ 6 ]. Early modeling efforts utilized the language of continuum mechanics to explain changes that occur on the surface (appositional growth) or throughout volumes of living tissues such as plants and bones [ 7 ]. Later, these approaches were supplemented with the finite element method to capture the spatially discontinuous aspects of differential growth [ 6 ]. One of the most notable growth theories is kinematic growth (KG) theory [ 6 , 7 ]. It posits that growth can be modeled through the evolution of the local growth rate tensor, which defines the magnitude and orientation of tissue expansion (or resorption). As neighboring subvolumes grow into each other, elastic stresses emerge to maintain subvolume boundaries. Rodriguez et al. extended KG by formulating the deformation gradient via the multiplicative split of the growth tensor and an elastic tensor [ 8 ]. They demonstrated that residual stresses in arterial tissue emerge due to incompatible growth between otherwise stress-free adjacent tissue subvolumes. Further, they modeled stress-dependent femur growth through a stress-sensitive growth tensor. KG has been popular due to its conceptual simplicity where changes in tissue volume and structure can be related to cell-level behaviors by making assumptions about how the growth tensor evolves in response to local state variables. Another popular approach was recently proposed for tissue adaptation based on constrained mixture (CM) theory [ 9 ]. The mass fractions of extracellular matrix constituents within tissue subvolumes (e.g., collagen, elastin, and glucosaminoglycans) or cell populations (modeled as concentrations) evolve over time. Tissue deformation results from changes in the mass of matrix components or changes in fluid pressure from osmotic effects, charged matrix constituents, or remodeling of tissue architecture (e.g., porosity) [ 10 ]. Mechanosensitivity may be introduced by assuming stress- or strain-sensitive production/conversion/degradation rates of matrix components or cell concentrations. The popularity of CM theory is only undercut by its relative complexity when compared to KG. CM theory requires tracking of the evolving mass fractions, introduces other constraints, increasing the degrees-of-freedom, computational expense, and effort to linearize finite element solutions. Kinematic growth and constrained mixture theories have each successfully been used to explore stress- and strain-driven tissue growth in contexts including morphogenesis, cardiovascular disease, cancer, and surgical planning [ 2 , 11 – 15 ]. A continuing challenge has been coupling mechanosensitive, biochemical cell signaling with tissue-level G&R. In the past, cell, and tissue growth were formulated in direct response to mechanical stimuli, neglecting the underlying mechanosensitive biochemical regulatory pathways and transport that govern cell and tissue growth [ 16 ]. New systems biology approaches are exploring growth in response to cell-level phenomena (see reviews by Refs. [ 5 ], [ 17 ], and [ 18 ]). Cell-scale and subcellular activities are regulated by their mechanical environment via mechanotransduction [ 19 ]. In response to mechanical stimuli, cells grow, migrate, differentiate, remodel the matrix, etc. Local and tissue-level mechanical stimuli cause cells to alter the microscale biochemical environment in a closed feedback loop to either restore mechanical homeostasis or bring the system into a new homeostatic state [ 20 – 23 ]. At the tissue-scale, cellular scale activities are often manifested as addition of mass (e.g., tissue thickening) and changes in material properties (e.g., stiffness and fiber organization). The integrated effects of these multiscale, multiphysical interactions have been captured by agent-based modeling on populations of individual cells or by modeling cell populations as continuous fields [ 11 , 24 – 27 ]. With more recent interest in systems modeling approaches [ 28 ], several challenging modeling problems arise, such as how spatial and temporal discretization can be resolved at the cell and tissue scales. Another open question concerns the choice between one- or two-way coupling of tissue mechanics and cell signaling. In one-way coupling, the cell-scale model drives G&R without any feedback. However, in two-way coupling, the cell-scale model is assumed to be mechanosensitive and changes in the tissue-scale mechanical environment as a consequence of G&R can feedback to affect the cell-scale model. The cell-scale systems biology model can be described with either systems of ordinary differential equations (ODEs) or PDEs, the latter accounting for spatial variation in species via processes such as diffusion. It is not clear how these modeling choices affect the multiscale G&R model. Most computational models of G&R coupled with systems biology have relied on commercial software [ 27 ], which is often difficult to customize for specific research needs. Commercial platforms have been used to model G&R, but they often require user-defined subroutines (such as abaqus ). Infrequently, custom in-house code may be provided, which is a less desirable option for the purposes we have outlined. Moreover, user subroutines for commercial software and in-house codes have no guarantee of forward compatibility. The use of commercial and in-house software limit innovation and reproducibility and creates barriers for researchers or educators who may not have the resources or training to access or modify these tools. The interdisciplinary nature of multiscale G&R highlights the clear need for open, transparent, and customizable computational tools to support collaborative research. Two notable open-source finite element software packages are common within the biomechanics community– febio (finite elements for biomechanics and biophysics [ 29 ]) 3 and fenics 4 [ 30 ]. febio is tailored for the biomechanics community, with relevant material constitutive models, boundary conditions, and physics built in. Users can define and run problems within febio 's friendly graphical user interface (GUI). febio is extendable through user-generated dynamically linked libraries built from C++ code [ 31 ]. Recently, a model repository was added that enables users to submit their published models. Currently, FEBio supports growth based on constrained mixture theory. A basic implementation of kinematic growth is provided, which users have previously extended through the plugin interface [ 32 ]. fenics is a general-purpose software for solving partial differential equations and can be imported as a module in python or built from C++ code. fenics takes advantage of just in time compilation (JIT) and autodifferentiation, allowing users to define problems without extensively deriving tangents used during linearization of the nonlinear finite element weak forms. On occasion, fenics have been used to simulate kinematic growth and remodeling [ 14 , 33 ]. Here, we develop open-source software for multiscale simulations of kinematic G&R of soft tissues driven by mechanosensitive signaling networks of chemical species described using the standard approach in systems biology. This effort is motivated by the fundamental need to integrate mechanosensitive systems biology with simulations of multiscale G&R. We view a mixture-level balance of the mass equation as central to coupling outputs of cell signaling to mathematical formulations for tissue growth. As a first step toward this goal, we focus on KG theory over CM theory because KG theory is more conceptually straightforward, is less computationally demanding, and CM theory is already sufficiently available within febio . We demonstrate the software capabilities through examples of (1) growth of a tissue graft with systems biology modeled by mass-action kinetics, and (2) aneurysm growth with systems biology driven by Hill-type kinetics. These problems are solely illustrative and detailed biological and physical aspects found in the literature are not considered here for simplicity. For each case, we compare one-way coupling (morphoelastic; i.e., biochemical signaling influences the growth formulation) and to-way coupling (biochemical signaling parameters are mechanosensitive and influence the growth formulation). The biochemical signaling for each example features position-dependent aspects to demonstrate growth heterogeneity. Additionally, the differences in ODE and PDE approaches for systems biology modeling are compared. With the goal of providing free, user-friendly, open-source software implementations to the biomechanics and mechanobiology communities, we implement and cross-verify our solutions in the open-source finite element software packages febio and fenics . 2 Methods We first review the kinematic growth theory within the context of solid tissue mechanics. Second, we define the governing equations for chemical evolution. This includes statements for mass-action chemical kinetics of molar concentrations, logic-based kinetics for normalized activity levels, and diffusive transport (for systems of PDEs). Third, we provide an overview of software implementations for febio and fenics . Finally, we describe the benchmark test problems. 2.1 Kinematic and Mechanics 2.1.1 Kinematics of Deformation. In continuum mechanics, constitutive models relate the energy in the form of stress that is required to deform a material from an initial to a deformed configuration. The motion χ maps position vectors X associated with the initial (material) configuration Ω 0 to position vectors x in the deformed (spatial) configuration Ω via x = χ ( X , t ) . Thus, the deformation gradient is a two-point tensor given by F = ∇ χ = ∂ x ∂ X ≡ Grad x (1) where Grad ≡ ∂ / ∂ X . Kinematic growth theory assumes that growth results from separate, successive deformations so that the deformation gradient can be multiplicatively decomposed as [ 8 , 34 ] F = F e F g (2) where the first deformation gradient ( F g ) occurs due to growth, and yields a fictitious “intermediate configuration” Ω * characterized by overlapping regions of newly grown material. The specific formulation of F g depends on a user-defined relationships between a stimulus and the rate of growth (detailed below). Adjacent material volumes may overlap in the intermediate configuration (i.e., the growth of one subregion is not constrained by the boundaries of its neighbors). The elastic deformation F e maps kinematic quantities from the intermediate configuration Ω * to the final deformed configuration Ω . This accounts for surface or body loads, as well as the stresses needed to resolve material overlap introduced in the intermediate configuration. Kinematic growth theories are appealing due to the ease of formulating the growth tensor as a function of any stimulus. The simplest formulations constrain F g to an orthonormal, symmetric positive-definite tensor of the form F g = ∑ i = 1 N ϑ i n i ⊗ n i , N ∈ { 1 , 2 , 3 } (3) where ϑ i are scalar growth stretches along the orthogonal directions n i [ 35 ]. In the absence of growth, ϑ i = 1 and F g = I where I is the identity tensor. Three types of growth can be produced from these constraints: line- ( N = 1 ), area- ( N = 2 ), and volume- ( N = 3 ) type growth. Other formulations may be explored at added computational and theoretical costs [ 34 , 36 ]. The specification of F g is assumed from measurements of tissue-level growth and the organization of individual tissue constituents (e.g., matrix fibers, cardiomyocytes, etc. [ 35 , 37 ]). In this work, the rate of growth ϑ ˙ is postulated to be given by ϑ ˙ = k ( ϑ ) ϕ ( • ) (4) where k ( ϑ ) is a function restricting unlimited growth, and ϕ ( • ) is the environmental function that scales the rate of growth in response to local conditions (e.g., stress, molar concentrations, etc..) [ 27 , 35 , 37 , 38 ]. The growth rate and environmental function reflect bioregulatory mechanisms that modulate growth in response to disruptions or deviations in state (e.g., mechanical stress). These mechanisms tend to return to homeostatic mechanical states but may also tend toward maladaptive states in response to disease. We implemented k ( ϑ ) as a bandpass function given by the sum of oppositely oriented sigmoidal functions k ( ϑ ) = ϑ min + ϑ max ( 1 1 + exp ( − ϑ − ϑ a γ ) + 1 1 + exp ( − ϑ + ϑ a γ ) ) (5) where ϑ min , ϑ max , ϑ a , and γ are user-defined parameters. For this study, the environmental function ϕ is related to the chemical kinetics of a biochemical species α with molar concentration (per unit volume in the reference configuration) denoted as c α . Specific formulations for k ( ϑ ) and ϕ are provided in the description of each test problem. The current growth stretch at each time is approximated through a Euler scheme as ϑ n = ϑ n − 1 + Δ t ϑ ˙ . 2.1.2 Mechanics of Kinematic Growth. To solve for the unknown deformation gradients F and F e , we first postulate F g which maps from the referential configuration to the intermediate configuration. The intermediate configuration is considered stress-free since growth is assumed to occur over a much longer time scale than subsequent elastic deformations and the daily energy needed for the thermodynamics of growth is negligible [ 39 , 40 ]. F e maps from the intermediate configuration to the deformed configuration and introduces the elastic stress tensor S e associated with the intermediate configuration Ω * . The second Piola–Kirchhoff stress tensor is related to the elastic stress by the pull-back S = J g F g − 1 S e F g − T where J g = det F g . The Cauchy stress follows via the push-forward operation [ 34 ]: σ = 1 J F S F T = 1 J e F e S e F e T (6) As a consequence of timescale differences, kinematic growth theories postulate the strain-energy function W e to capture the strain energy required to move from the intermediate configuration Ω * to the deformed configuration Ω . To evaluate W e , we define the elastic and growth right Cauchy-Green deformation tensors as C e = F e T F e and C g = F g T F g , respectively. Similarly, the change in volume J = det F (i.e., the Jacobian) can be split into elastic and growth contributions J e and J g via J = ( det F e ) ( det F g ) = J e J g . In this study, the strain energy W e was given by the unconstrained, coupled hyperelastic neo-Hookean constitutive model W e ( C e ) = c 1 ( I 1 e − 3 ) − 2 c 1 ln J e + K 2 ( ln J e ) 2 (7) where c 1 is a material constant, I 1 e = tr ( C e ) is the first invariant of the right Cauchy–Green strain tensor, and K is a material constant similar to the bulk modulus. The second Piola–Kirchhoff stress in the intermediate configuration is computed from the elastic strain energy via S e = 2 ∂ W e ∂ C e (8) The Cauchy stress follows from Eq. (6) . 2.2 Formulation of Chemical Kinetics in Systems Biology. We base our formulation of chemical kinetics on the works of Refs. [ 41 ] and [ 42 ]. Soft biological tissues are typically hydrated materials containing a solid phase and an interstitial fluid phase. Chemical species may move with either the solid or the fluid phase. A comprehensive framework would need to account for the evolution of these phases as well as their volume fractions, orientations, and the diffusive drag between phases (i.e., the hydraulic pressure of the interstitial fluid). For the present study, we relax these assumptions since their contributions are minimized by the difference in time-scales between mechanical loading, tissue growth, and chemical kinetics. We assume that solid tissue deformation does not directly affect chemical kinetics, we do not explicitly model the interstitial fluid, and we assume ideal behavior of chemical species. More information about these assumptions is provided in Sec. A.1 . 2.2.1 Mass-Action Kinetics. A generic mass-action chemical reaction is given by ∑ ν R α c α ↔ ∑ ν P α c α (9) where the molar concentration of a chemical species α is given by c α , ν R α is the stoichiometric coefficient of species α as a reactant, and ν P α is the stoichiometric coefficient of species α as a product [ 41 , 42 ]. The net stoichiometric coefficient for species α is thus ν α = ν P α − ν R α . The molar supply c ^ i α is related to the production rate ζ ^ i of chemical reaction i according to c ^ i α = ν i α ζ ^ i (10) The associated mass-action molar production rates are given as ζ ^ F = k F ∏ α ( c α ) ν R α , ζ ^ R = k R ∏ α ( c α ) ν P α , ζ ^ = ζ ^ F − ζ ^ R = ζ ^ F [ 1 − k R k F ∏ α ( c α ) ν α ] (11) for forward, reverse, and reversible reactions, respectively. Implicitly available concentration supplies ( c ^ q α ) and sinks ( c ^ d α ) are modeled with the mass-action kinetics given when ν i α = 1 c ^ q α = k q , c ^ d α = − k d c α (12) 2.2.2 Logic-Based Reaction Kinetics. The kinetics of complex biochemical signaling networks may be simplified via logic-based governing equations in the form of normalized Hill equations as demonstrated by Ref. [ 12 ]. Here, concentrations or activity levels of biochemical species are normalized in the range c α ∈ [ 0 , 1 ] . The concentration/activity levels evolve through activation or inhibition by other signaling members which are governed by a Hill function f ∈ [ 0 , 1 ] . The molar production rate is found by scaling the Hill function f with the weight k H c ^ i α = k H f i α (13) where i indicates a unique reaction. First, we consider the production rate c ^ → β α which represents the production rate of a species α in response to activation ( → ) by a secondary species β . The normalized Hill activation function is given by f → β α = B ( c β ) n K n + ( c β ) n (14) where B = ( EC 50 n − 1 − 1 ) / ( 2 EC 50 n − 1 − 1 ) , K = ( B − 1 ) ( 1 / n ) , n is the Hill coefficient, and EC 50 is the value of c β at which half of the maximal activation of α occurs. The inhibition ( ¬ ) of species productions is determined via Boolean logic. For example, consider the normalized production rate of a species α that can be inhibited by a tertiary species γ c ^ ¬ γ α = k H ( 1 − f → γ α ) = k H f ¬ γ α (15) where ¬ is the Boolean NOT operator. The Boolean AND operator ( ∧ ) allows us to consider the simultaneous activation of α by β and inhibition of α by γ c ^ → β ∧ ¬ γ α = k H ( f → β α ( 1 − f ¬ γ α ) ) (16) 2.2.3 Chemical Transport. In this study, chemical transport occurs solely due to diffusion, i.e., j α = − D α ∇ c α = − d α · grad c α , ∂ c α ∂ t = D α ∇ 2 c α = grad c α · d α · grad c α (17) where j α is the molar flux, d α = D α I is the isotropic diffusivity tensor, and grad ≡ ∂ / ∂ x . The general expression for the total molar change is given by the sum of diffusive and reactive terms for each chemical species α and each chemical reaction i , i.e., d c α d t = D α ∇ 2 c α + ∑ i c ^ i (18) 2.3 Software Implementation and Solution Strategies 2.3.1 Code Implementation. The G&R models are implemented in the open-source finite element software packages febio [ 29 ] and fenics [ 30 ]. febio is available as an executable application alongside febiostudio —a graphical user environment that can be used for all steps in the simulation pipeline including model generation, solution, data analysis, and visualization. fenics is a general tool for solving PDEs governed by weak forms defined in symbolic expressions using the fenics Form Language. The input files for the test problems in this study are available on GitHub at the following links. 5 , 6 Mechanosensitive KG is implemented in febio version 4.7 as a custom plugin (dynamically linked library) that can be imported into febio [ 31 ]. The plugin extends the standard febio solver modules to introduce Hill-type kinetics and stress-dependent KG and reaction kinetics. The plugin also relaxes some default assumptions (see Sec. A.4 for full details). Model geometry, meshing, physics, and coefficients are specified through an xml-formatted text file. Users can generate models and assign coefficients through the FEBioStudio graphical interface or manually edit input files. Model results can be visualized in FEBioStudio or exported to VTK format for visualization in paraview . Mechanosensitive KG is implemented in fenics ( fenics 2019.1.0, Python 3.6). A python script establishes the mechanical and chemical equations and a main subroutine to step through the solution and save model outputs, and read model inputs (the geometry, mesh, and coefficients). Model results are exported to xdmf/h5 formats for visualization in paraview . 2.3.2 Solution Strategies. Although the same statement of virtual work applies to models in febio and fenics , the implementation strategies used to linearize and discretize the solutions differ as summarized in Fig. 1 . In febio , we employed a single solver to a monolithic problem where the linearization and discretization are explicitly derived and hard-coded. In fenics , we employed separate solvers to evaluate nodal displacements u and nodal chemical concentrations { c α , … , c γ } T = c , and the linearization and discretization were achieved through autodifferentiation functions. Fig. 1. Open in a new tab Algorithms used in febio and fenics . The solution in febio is monolithic (i.e., nodal displacements u and nodal concentrations c are solved from a single virtual work equation). The solution in fenics first evaluates virtual work for nodal concentrations c and then solves for nodal displacements u . One-way coupling and two-way coupling are applied after by updating s ( σ ) after solving u and c . For both cases, we explore one-way and two-way coupling of the systems biology and mechanics (Fig. 1 ). The one-way coupling mechanism assumes that growth is solely driven by biochemical signaling where reaction rates and prescribed concentrations are independent of mechanical forces. The to-way coupling mechanism assumes that reaction rates or prescribed concentrations are reciprocally sensitive to stress (enforced by the scaling functions s ( σ ) defined in Secs. 2.4.3 and 2.5.3 ). 2.3.3 Finite Element Specification. Tetrahedral elements are used for all meshes. In fenics , we are limited to four node tetrahedrons (tet-4, on-vertex nodes) equipped with quadratic interpolation and straight element edges. In febio , 10-node tetrahedrons (tet-10, on-vertex, and midedge nodes) are equipped with quadratic interpolation and curved element edges. 2.4 Test Problem 1: Axial Growth of a Tissue Graft 2.4.1 Problem Overview. The first test problem demonstrates tissue-level differential growth of an idealized cubic tissue graft over 10 years (Fig. 2( a ) ). Growth occurs downstream of nutrients produced by a synthesis reaction and implicitly generated supplies as modeled by spatially varying mass-action kinetics. We model mechanosensitive growth by scaling implicit production rates dependent on local stress. Finally, we demonstrate the effects of diffusive transport on growth for the one-way cases. Fig. 2. Open in a new tab Test problem descriptions. ( a ) A cubic tissue graft undergoes growth in the vertical (Z) direction. A 1 4 symmetry mesh is used. Biochemical species A has a spatially varying production rate given by c ^ A ( X , s ) . The growth rate is proportional to biochemical species C. ( b ) A vessel experiences area-type growth (X, Y) in response to matrix metalloproteinases (MMP) downstream of TGF β . TGF β is spatially distributed according to the boundary condition c max TGF β ( X , s ) . The signaling network is governed by logic-based kinetics (arrow: activation, blunted arrow: inhibition). The growth rate is proportional to the rate of MMP activity, which is inversely proportional to collagen concentration. 2.4.2 Graft Geometry and Displacement Boundary Conditions. The graft is idealized as a unit cube ( 1 mm 3 , Fig. 2( a ) ). A 1 4 symmetry geometry is employed to model only one quadrant of the graft (split along the X- and Y -axes). Symmetry boundary conditions are applied by constraining the displacement of surfaces coinciding with the x -axis to the XZ plane and those of the y -axis to the YZ plane. The bottom of the graft is fixed. 2.4.3 Graft Systems Biology. A system of differential equations with arbitrary chemical species A, B, and C is modeled with chemical reactions, implicit concentration supplies, and implicit species degradation governed by A + B → k F C , Mass action reaction (19) A , B , C → k d α [ ] , Implicit degradation (20) [ ] → k q α A , B , C . Implicit supply (21) To illustrate differential growth (e.g., due to a spatial distribution of solid-bound enzymes, cells, or other upstream stimulus), the production rate of one species for each case is specified as a function of position k q ( X ) (one-way coupling) or a function of position and stress k q ( X , σ ) (two-way coupling). The positional dependence is introduced through the function q graft ( X ) = exp ( − 10 3 ( | X | + | Y | ) 2 ) (22) where X and Y are coordinates in the reference configuration (Fig. 2( a ) ). For two-way coupling, the Cauchy stress is calculated after each step, and the term k q ( X , σ ) is updated based on this stress. These position-dependent production rates may represent a heterogeneous quality of the tissue (e.g., distribution of solid-bound enzymes, cells, or some upstream stimulus). The production rates are formulated as a function of the initial two-dimensional position given by k q ( X ) = q graft ( X ) , one-way coupling (23) k q ( X , σ ) = q graft ( X ) s ( σ ) , two-way coupling (24) where s ( σ ) is a sigmoid function dependent on the trace of Cauchy stress tr ( σ ) given by s ( σ ) = a 0 + a amp 1 + exp ( − tr ( σ ) − s 0 b ) (25) with a 0 = 1.0 , a amp = 10.0 , b = 0.05 , and s 0 = 0.3 . The parameters for the bandpass function 5 are set to ϑ min = 0.0 , ϑ max = 1.0 , ϑ a = − 3.0 , and γ = 0.3 . The initial concentration for each species is 0 mM. The resulting governing equations are summarized as d c A d t = D ∇ 2 c A − k F c A c B + k q A − k d A c A (26) d c B d t = D ∇ 2 c B − k F c A c B + k q B − k d B c B (27) d c C d t = D ∇ 2 c C + k F c A c B + k q C − k d C c C (28) Four different cases are presented for test problem 1 to illustrate the influence of coupling schemes, diffusivity, and spatial heterogeneity inherent to biological tissues. The parameters for these cases are listed in Table 1 . Different diffusivities D ∈ { 0 , 1 × 10 − 10 , 1 × 10 − 9 , 1 × 10 − 8 } mm 2 / year are evaluated for the one-way coupled systems. Table 1. Kinetic parameters used for the tissue graft problem cases k F k q A k q B k q C k d A k d B k d C Case Coupling ( mM · year ) − 1 ( mM / year ) ( mM / year ) ( mM / year ) ( 1 / year ) ( 1 / year ) ( 1 / year ) 1 One-way 1 k q ( X ) 0.25 0 1 0.25 0.25 2 Two-way 1 k q ( X , σ ) 0.25 0 1 0.25 0.25 3 One-way 1 0.25 1 k q ( X ) 0.25 0.25 1 4 Two-way 1 0.25 1 k q ( X , σ ) 0.25 0.25 1 Open in a new tab 2.4.4 Kinematics of Graft Growth. The growth tensor is constrained to line-type growth in the vertical direction [ 34 , 43 ] F g = I + [ ϑ − 1 ] n Z ⊗ n Z , (29) where ϑ is obtained using Eq. (4) with the environmental function dependent on species C defined as ϕ ( c ˙ C ) = c ˙ C (30) with c ˙ C approximated by a backward Euler scheme. 2.5 Test Problem: Aneurysm Growth 2.5.1 Problem Overview. The second test problem demonstrates tissue-level aneurysm growth in response to transforming growth factor beta (TGF β ) signaling, inspired by examples in the literature [ 12 , 44 , 45 ]. Tissue growth occurs downstream of locally elevated TGF β over 8 years of growth. The signaling network dynamics are governed by logic-based kinetics (Fig. 2( b ) ). We model mechanosensitive growth by scaling TGF β activity in response to local stresses. Finally, we demonstrate the effects of diffusive transport on growth for the one-way cases. 2.5.2 Aneurysm Geometry and Displacement Boundary Conditions. An initially healthy vessel is represented by a 6 mm long tube with 0.5 mm inner and 0.65 mm outer diameters. A 1 8 symmetry geometry (along X -, Y -, Z -planes) is represented with one octant of a tube with symmetry planes in the XY- , XZ- , and YZ -planes (Fig. 2( b ) ). Symmetry boundary conditions are similarly imposed by fixing the displacement of surfaces coinciding with along the X -axis to the XZ plane, those of the Y -axis to the YZ plane, and those of the Z -axis to the XY plane. Diffusive signal member fluxes are zero along the normals of symmetry faces. The vessel is pressurized to 120 mmHg at the beginning of each simulation and is constantly applied over the course of eight years of growth. 2.5.3 Aneurysm Systems Biology. Aneurysm growth is a complicated process involving signaling between numerous biochemical pathways including TGF β . We model a simple TGF β -driven signaling network graphically described in Fig. 2( b ) based on Ref. [ 12 ]. The primary species TGF β binds to cell surface receptors, which activates intracellular mitogen-activated protein kinase (MAPK) or Smad. MAPK and Smad translocate to the nucleus where MAPK activates expression of genes for matrix metalloproteinases (MMP) and Smad activates the tissue inhibitors of metalloproteinases (TIMPs). MMPs and TIMPs affect the degradation of the ECM and may release matrix-bound growth factors including TGF β . The activity of all species except for TGF β is modeled to decay due to various biophysical means (species degradation, competition, and receptor recycling, etc.). Due to the complexity and uncertainty of biochemical signaling and experimental measurements, logic-based kinetics are used to describe the governing equations. In a continuum, the activity of signaling members may be considered a normalized concentration reflecting the average activity level of local cells. For this problem, we assume that TGF β activity is elevated near the center of a vessel (due to stress, disease, etc.). MAPK and Smad are intracellular members and do not diffuse in our model. MMP and TIMP are soluble products and may diffuse. The position-dependent TGF β activity for this problem (Fig. 2( b ) ) is given by q vessel ( X ) = exp ( − | Z | 3 3 ) (31) For the two-way coupling method, we evaluate the effect of mechanosensitive TGF β activity by calculating the Cauchy stress after each step and prescribing TGF β activity ( c TGF β ) based on this stress. The governing equations are detailed as follows: c TGF β = { c max TGF β q vessel ( X ) one-way coupling c max TGF β q vessel ( X ) s ( σ ) two-way coupling (32) d c MAPK d t = k H f → TGF β MAPK − k d c MAPK (33) d c Smad d t = k H f → TGF β Smad − k d c Smad (34) d c MMP d t = D ∇ 2 c MMP + k H f → TGF β ∧ ¬ TIMP MMP − k d c MMP (35) d c TIMP d t = D ∇ 2 c TIMP + k H f → Smad TIMP − k d c TIMP (36) where n = 1.2 , EC 50 = 0.5 , a 0 = 1.0 , a amp = 1.0 , b = 0.2 , and s 0 = 0.05 . The degradation rate is k d = 0.25 year − 1 for all species except the solid-bound TGF β . One-way and two-way coupling is introduced through the prescribed initial condition for c TGF β . The one-way coupling approach defines c TGF β solely by the position along the axis of the vessel (coinciding with the Z direction; Z = 0 is the center of the vessel). The two-way coupling approach additionally defines c TGF β as a function of the mechanosensitive sigmoidal function s ( σ ) . The same parameters for s ( σ ) and k ( ϑ ) were used as in test problem 1. The initial concentration of all other species is zero. Three cases are presented to illustrate the influence of coupling schemes, diffusivity, and spatial heterogeneity inherent to aneurysm growth. Case 3 also investigates the effect of stimulus intensity by doubling the prescribed initial concentration c TGF β . The parameters for these cases are outlined in Table 2 . Diffusivity values D ∈ { 0 , 1 × 10 − 10 , 1 × 10 − 9 , 1 × 10 − 8 , 1 × 10 − 7 , 1 × 10 − 5 } mm 2 / year are evaluated for MMP and TIMP for case 1. Table 2. Kinetic parameters and prescribed TGF β activity for the aneurysm problems Case Coupling k H ( 1 / year ) c TGF β ( X ) c max TGF β 1 One-way 1 c TGF β ( X ) 0.5 2 Two-way 1 c TGF β ( X ) 0.5 3 One-way 1 c TGF β ( X ) 1 Open in a new tab 2.5.4 Kinematics of Aneurysm Growth. Area-type growth is constrained along the radial direction of the vessel (i.e., normal to the axial direction n Z of the vessel) [ 43 ] F g = ϑ I + [ 1 − ϑ ] n Z ⊗ n Z . (37) MMP simultaneously degrades collagen and releases matrix-bound growth factors which contribute to delayed collagen deposition. This results in both decreases and increases in total collagen, yet MMP concentration is overall correlated with tissue growth [ 46 ]. This discrepancy is explained by either a decrease in the collagen volume fraction (resulting in tissue-level volumetric increases due to differences in density between individual tissue constituents and increased permeability of degraded tissues) or by the delayed addition of new collagen mass due to downstream effects of releasing matrix-bound growth factors [ 47 – 49 ]. To keep this illustrative model simple, we focus on the correlation between MMP concentration and aneurysm growth. The downstream processes are internally captured through ϕ ( c ˙ MMP ) rather than explicitly tracking evolving collagen fractions ϕ ( c ˙ MMP ) = 6.5 d c MMP d t (38) where c ˙ MMP is approximated by a backward Euler scheme. 2.6 Cross-Verification and Mesh Discretization Study. Cross-verification is performed by comparisons between febio and fenics results for fields associated with growth (relative volume), mechanics (Lagrange strain magnitude, trace of the Cauchy stress tensor), and reactive species levels. Test problems are qualitatively and quantitatively compared with paraview 5.11 (Kitware, Clifton Park, NY). Solutions from febio tet-10 meshes are projected onto fenics tet-4 meshes. Data fields from both software are appended to the new mesh and the absolute difference is calculated for relevant fields. The deformed fields are visualized with the average displacement. Due to the highly skewed nature of difference distributions, a blue-grey-red heat map is modified to correspond to the 5th, 50th, and 95th quantiles of the absolute differences. Additional quantitative comparison is achieved through a mesh convergence study. The convergence of each software in response to one-way or two-way coupling is evaluated by mesh refinement (number of elements) and biasing (i.e., linearly varying the spatial interval between nodes to produce smaller or larger elements without affecting the total number of elements). Detailed mesh convergence results are provided in Supplemental S1 available in the Supplemental Materials on the ASME Digital Collection. Test problem 1 (tissue graft) is discretized with 120, 324, 672, 1,200, 9,600, or 32,400 elements (see Supplemental Figure S1 available in the Supplemental Materials ). We also investigate biasing elements to be smaller in regions of high growth in the center (0,0,1) and at the corner (0.5, 0.5, 1). There are a total of 12 meshes with six levels of element discretization and two levels of biasing (unbiased and biased). Convergence of the relative volume, trace of stress, and concentration c C are compared for each mesh and software at the top surface of the graft. Test problem 2 (aneurysm) is discretized with 450, 1,920, 5,400, 8,100, 24,000, or 52,500 elements (see Supplemental Figure S6 available in the Supplemental Materials on the ASME Digital Collection). We also consider biasing elements to be smaller in regions of high growth near the center of the Z axis (plane normal to 0,0,1) and toward the intimal and abluminal walls ( r = 0.5, 0.65). This results in 12 meshes similar to the first test problem. Convergence of the relative volume, trace of stress, and c MMP are compared for each mesh and software at the luminal surface near r = 0.5, Z = 0 . 3 Results 3.1 Test Problem 1: Tissue Graft. Figures 3 and 4 show growth patterns for graft cases 1–4 at different times with contour plots of the volume ratio (i.e., J = V V 0 ) (panels A). Cases 1 and 2 appear to undergo constant growth over the course of 10 years. In contrast, cases 3 and 4 grow rapidly within the first 3 years but grow minimally afterward. Endpoint comparisons between febio and fenics are displayed for the volume ratio, trace of stress ( tr ( σ ) ), and the concentration of the morphogen species [ C ] (panels B-D). An exposed view is also included to visualize internal distributions. Fig. 3. Open in a new tab Results for test problem 1 (cubic tissue graft) cases 1 (one-way coupling) & 2 (two-way coupling) with no diffusivity (ODE chemical system). Growth results downstream of biochemical species C. The supply rate of species A is position-dependent. ( a ) Tissue growth (volume ratio) occurs over a period of 10 years ( febio results shown). ( b – d ) Comparison between febio and fenics of volume ratio, trace of Cauchy stress, and concentration of chemical species C. A surface view and a 3 4 symmetry exposed view is presented to visualize the external and internal growth. Simulation results are qualitatively nearly identical. The absolute differences in distributions between febio and fenics data fields are quantified (right column of ( b – d )) with color maps that vary to indicate quantiles ( q ) that correspond to 5% (blue), 50% (i.e., the median and gray), and 95% (red). ( a and b ) Growth is more homogenous in the two-way case than the one-way case. The greatest growth is observed in the center of each case. ( c ) The two-way coupling facilitates negative feedback, reducing stress throughout the graft. ( d ) Production of species C increased at the center and top of the graft due to two-way coupling. Quantified differences are generally orders of magnitude below the range for the corresponding data fields, with no clear differences appearing. (Color version online.) Fig. 4. Open in a new tab Results for test problem 1 (cubic tissue graft) cases 3 (one-way coupling) & 4 (two-way coupling) with no diffusivity (ODE chemical system). Growth results downstream of biochemical species C. The supply rate of species C is position-dependent. ( a )Tissue growth (volume ratio) occurs over a period of 10 years ( febio results shown). ( b – d ) Comparison between febio and fenics of volume ratio, trace of Cauchy stress, and concentration of chemical species C. A surface view and a 3 4 symmetry exposed view are presented to visualize the external and internal growth. Simulation results again are qualitatively nearly identical. The absolute differences in distributions between febio and fenics data fields are quantified (right column of ( b – d )) with color maps that vary to indicate quantiles ( q ) that correspond to 5% (blue), 50% (i.e., the median, gray), and 95% (red). ( a and b ) Growth is again more homogenous in the two-way case than the one-way case. The greatest growth is observed in the center of each case. ( c ) The two-way coupling facilitates negative feedback, reducing stress throughout the graft. ( d ) Production of species C increased at the center and top of the graft due to two-way coupling. Growth for cases 3–4 is greater than for cases 1–2, with the primary difference being which species has a position-dependent supply given by k q ( X ) or k q ( X , σ ) . Differences are orders of magnitude below the corresponding data fields. For case 4, larger differences in tr ( σ ) and [C] appear along sharp gradients near the top surface, indicating that febio and fenics predict slightly different gradient boundaries. (Color version online.) The two-way coupling schemes in cases 2 and 4 leads to modifications in patterns of growth when compared to cases 1 and 3. Notably, the two-way coupling models illustrate negative feedback where stress amplifies tissue expansion in regions of high stress. This results in more homogeneous distributions of growth, ultimately reducing stress and strain gradients. The influence of chemical diffusivity for cases 1 and 3 (one-way) is illustrated in Figs. 5 and 6 . Tissue growth, stresses, and biochemical species concentrations become more homogeneous with increasing diffusivity D . At high enough diffusivities, growth becomes nearly uniform and delocalized from the position-dependent production functions k q ( X ) A and k q ( X ) C . The chemical kinetics of species C (which drives growth for each case) over time is summarized by the maximum and minimum concentrations associated with each model and diffusivity value in Fig. 7 . Chemical equilibrium is reached shortly after 10 years for cases 1 and 2 but around 4 years for cases 3 and 4. Equilibrium took longer to achieve for cases 1 & 2 since the position-dependent source k q ( X ) A was upstream of the morphogenic species C. With increasing diffusivity, the maximum and minimum values begin to converge toward the average concentration. Fig. 5. Open in a new tab Selected results for test problem 1 (cubic tissue graft) case 1 (one-way coupling) with diffusivity (PDE). Growth results downstream of biochemical species C. The supply rate of species A is position-dependent. ( a ) Tissue growth (volume ratio) occurs over a period of 10 years ( febio results shown). ( b – d ) Comparison between febio and fenics of volume ratio, trace of Cauchy stress, and concentration of chemical species C. A surface view and a 3 4 symmetry exposed view are presented to visualize the external and internal growth. Simulation results are qualitatively nearly identical. The absolute differences in distributions between febio and fenics data fields are quantified (right column of B-D) with color maps that vary to indicate quantiles ( q ) that correspond to 5% (blue), 50% (i.e., the median, gray), and 95% (red). Growth (volume ratio), stress, and species C profiles become more spatially homogenous with increasing diffusivity. Differences are orders of magnitude below the corresponding data fields and decrease with moderate diffusivity. (Color version online.) Fig. 6. Open in a new tab Selected results for test problem 3 (cubic tissue graft) case 1 (one-way coupling) with diffusivity (PDE). Growth results downstream of biochemical species C. The supply rate of species C includes a position-dependent component. ( a ) Tissue growth (volume ratio) occurs over a period of 10 years ( febio results shown). ( b – d ) Comparison between febio and fenics of volume ratio, trace of Cauchy stress, and concentration of chemical species C. A surface view and a 3 4 symmetry exposed view are presented to visualize the external and internal growth. Simulation results are qualitatively nearly identical. The absolute differences in distributions between febio and fenics data fields are quantified (right column of ( b – d )) with color maps that vary to indicate quantiles ( q ) that correspond to 5% (blue), 50% (i.e., the median, gray), and 95% (red). Growth (volume ratio), stress, and species C profiles become more spatially homogenous with increasing diffusivity. Differences are orders of magnitude below the corresponding data fields and decrease with moderate diffusivity. (Color version online.) Fig. 7. Open in a new tab Maximum and minimum concentration of species [C] is plotted during the time for test problem 1 (cubic tissue graft) cases 1–4. Production of [A] is position dependent in cases 1–2. Production of [C] includes a position dependent term for cases 3–4. Diffusivity reduces heterogeneity in the distribution of each chemical species as well as downstream growth. Case 4 exhibits an order of magnitude increase in [C] compared to case 3. Growth in case 4 does not reflect this increase due to the bandpass function. Results from febio and fenics are nearly identical for each case. Quantitative analysis indicates that the solutions only differ for two-way coupled cases near sharp gradients in the solution (Figs. 3 and 4 ). Solutions for two-way coupled problems (cases 2 and 4) were more sensitive to mesh refinement and biasing elements toward regions of high growth (see Supplemental Figure S2–S5 available in the Supplemental Materials on the ASME Digital Collection). 3.2 Test Problem 2: Aneurysm. Figure 8 shows growth patterns for the aneurysm cases 1–3 at different times with contour plots of the volume ratio (i.e., J = V V 0 ) shown in panel A. An aneurysm forms for each case within 2 years, with case 3 predicting a larger aneurysm than cases 1 and 2. Aneurysms nearly triple in volume by year 8, with a double aneurysm emerging for case 3 before year 5. Endpoint comparisons between febio and fenics are displayed for the volume ratio, trace of stress ( tr ( σ ) ), and the normalized activity [MMP] of matrix metalloproteinases (panels B-D). An exposed view of the lumen is also included to visualize internal distributions after 8 years of growth. Fig. 8. Open in a new tab Results for test problem 2 (aneurysm) cases 1 (one-way coupling), 2 (two-way coupling), and 3 (one-way coupling, doubled TGF β activity) with no diffusivity (ODE chemical system). Growth results due to MMP activity downstream of TGF β . The supply rate of TGF β is position-dependent. ( a ) Aneurysm growth (volume ratio) occurs over a period of 8 years ( febio results shown). ( b – d ) Comparison between febio and fenics of volume ratio, trace of Cauchy stress, and MMP activity. Views of the abluminal and luminal walls are presented to visualize external and internal growth, respectively. Simulation results are qualitatively nearly identical aside from case 2 where the differences are limited to the abluminal surface. The absolute differences in distributions between febio and fenics data fields are quantified (right column of ( b – d )) with color maps that vary to indicate quantiles ( q ) that correspond to 5% (blue), 50% (i.e., the median, gray), and 95% (red). ( a and b ) Growth (volume ratio) changes only slightly on the abluminal wall due to two-way coupling. Doubling the prescribed TGF β in case 3 leads to a double aneurysm. ( c ) The max stress increases for case 2. The max and min stress become more extreme for case 3, indicating more tension and compression. ( d ) MMP activity redistributes in response to two-way coupling (case 2) and doubling the prescribed TGF β activity (case 3). Regions of MMP activity around 0.25 correlate to regions with the greatest growth. Differences are orders of magnitude below the corresponding data fields. For case 2, larger differences in tr ( σ ) and [C] appear along sharp gradients near the center of the abluminal surface. (Color version online.) Differences between the one-way and two-way coupling schemes in cases 1 and 2 are less visibly pronounced than in the graft test problem. In case 3 (one-way, doubled TGF β activity), a single, large aneurysm rapidly forms by 2 years. The aneurysm begins to split around 5 years and a double coexisting aneurysm forms by 8 years. The average measures for tr ( σ ) increased compared to cases 1 and 2. However, the maximum values of tr ( σ ) for one-way coupling (cases 1 and 3) are below those for two-way coupling (case 2), indicating more localized extreme tension due to two-way coupling. The influence of chemical diffusivity for case 1 (one-way) is illustrated in Figs. 9 and 10 . Tissue growth, stresses, and biochemical species concentrations become more uniform with increasing diffusivity D . At higher diffusivities, the vessel thickens rather than forming an aneurysm. Fig. 9. Open in a new tab Results for test problem 2 (aneurysm) case 1 (one-way coupling) with diffusive MMP and TIMP (PDE chemical system). Growth results due to MMP activity downstream of TGF β . The supply rate of TGF β is position-dependent. ( a ) Aneurysm growth (volume ratio) occurs over a period of 8 years ( febio results shown). ( b – d ) Comparison between febio and fenics of volume ratio, trace of Cauchy stress, and MMP activity. Views of the abluminal and luminal walls are presented to visualize external and internal growth, respectively. Simulation results are qualitatively nearly identical. The absolute differences in distributions between febio and fenics data fields are quantified (right column of ( b – d ) with color maps that vary to indicate quantiles ( q ) that correspond to 5% (blue), 50% (i.e., the median, gray), and 95% (red). Growth (volume ratio), stress, and MMP activity profiles become more spatially homogeneous with increasing diffusivity. With increased diffusivity, TGF β leads to vessel dilation/thickening rather than aneurysm growth. Differences are orders of magnitude below the corresponding data fields and decrease with moderate diffusivity. (Color version online.) Fig. 10. Open in a new tab Arterial thickness as a function of position along the Z axis (aligned with the axial direction) is plotted from the center to one end of the vessel after 8 years of growth for cases 1–3 (left, 1 2 geometry). The initial thickness was 0.15 mm. Variations in arterial thickness after 8 years of growth due to species diffusivity (right). High diffusivity leads to uniformly distributed growth. We examine the interesting emergence of the double aneurysm in case 3 by calculating the steady-state response of MMP to TGF β activity as done by Ref. [ 12 ]. The end-time profiles of TGF β and MMP are presented alongside the equilibrium sensitivity of MMP to TGF β in Fig. 11 . MMP has a nonmonotonic relationship with TGF β , and MMP activity peaks when TGF β activity is roughly half maximal. Thus, doubling the TGF β concentration in case 3 leads to the spatial delocalization of MMP activity and growth from TGF β . Fig. 11. Open in a new tab Doubling the TGF β concentration in case 3 (one-way, ODE) led to a double aneurysm. This occurs due to the nonmonotonic relationship between TGF β and MMP due to TIMP inhibition of MMP [ 24 ]. At equilibrium, MMP activity peaks at regions where TGF β ≈ 0.5 (as indicated by arrows). Results from febio and fenics are nearly identical for each case, with differences only occurring for case 2 (two-way) in regions with sharp gradients along the radial (transmural) direction (Figs. 8 and 9 ). Simulation results are relatively insensitive to mesh refinement for case 1 and case 3 only displays sensitivity with respect to the trace of stress (see Supplemental Figures S7–S9 available in the Supplemental Materials ). The two-way case (case 2) displays low levels of mesh sensitivity. Biasing the mesh toward the center or the lumen and abluminal wall had a minimal effect on solution values. fenics models failed iterative convergence for cases 1–3 using the coarsest mesh (450 elements) and case 2 for the second coarsest mesh (1,920 elements). The coarsest meshes only discretize the vessel wall with 1–2 elements, which in combination with tet-4 element formulations poorly capture transmural gradients. 4 Discussion Advances in continuum biomechanics modeling have enabled tremendous progress in computational modeling of G&R [ 43 ]. However, mechanistic modeling of biological growth is still an emerging field facing numerous hurdles. First, parameterization and quantitative characterization of cell signaling pathways are challenging, particularly once more detailed signaling networks are considered [ 50 ]. Second, the nature of G&R is a multiscale process [ 43 , 51 ], which necessitates that coupling methodologies accommodate different spatial and temporal scales. This requires better integration of experimental correlations between tissue-level growth and cell-level signaling pathways. Third, systems biology often models signaling pathways via normalized molar concentrations or activity levels while continuum mechanics tracks the real mass (density) of each constituent. These approaches help constrain chemical evolution and are useful when experimental measures for concentration and reaction kinetics are unknown or difficult to establish. Fourth, cell-level behavior (mechanosensitivity, signaling pathways, mutant phenotypes, etc.) is stochastic and some variation in cell-level behaviors is lost. Continuum approaches model the average tissue behavior and does not explicitly model cell growth and migration. Finally, the effective chemical concentration that is sensed at the cellular level may greatly differ from the concentration as measured per tissue-scale volume. This discrepancy emerges due to structural simplifications (e.g., tissue porosity), changes in the tissue volume due to elastic deformation, and osmotic effects. While these may be modeled within the framework of mixture theory, doing so introduces additional field variables (e.g., fluid pressure), greatly complicating linearization and discretization of the governing equations as well as computation time and resource management. In this study, we present a multiscale G&R model coupled with systems biology to study the interaction between cell-scale signaling events and tissue-scale growth. We adopt the growth rate Eq. (4) from Ref. [ 35 ] to drive growth based on chemical species rather than directly based on mechanical stress. This contrasts with most other kinematic growth modeling frameworks, where the growth tensor F g is historically related to the stress or strain [ 39 , 52 – 54 ] with fewer KG studies relating F g to chemically driven systems biology [ 27 , 49 , 55 ]. Reactive kinetics are modeled by mass-action or logic-based Hill-type kinetics, the latter of which has more recently been adapted for finite element modeling of G&R [ 2 , 12 , 32 , 56 ]. We model mechanosensitive homeostatic and maladaptive regulatory mechanisms through scaling of reaction kinetics by the stress-sensitive control function (Eq. (25) ). The illustrative examples presented here demonstrate potential roles of mechanosensitivity in growth processes. For the tissue graft problem, two-way coupling redistributes morphogenic chemical species and downstream tissue growth. Compared to one-way coupling, two-way coupling results in lower magnitude measures of stress and stress gradients characteristic of homeostatic negative feedback loops (Figs. 3( c ) and 4( c ) ). The effects of two-way coupling are more subtle for the aneurysm problem, where its effects are restricted to the abluminal surface. This may be due to implicit constraints on the chemical evolution due to normalized, logic-based chemical kinetics and inhibition of MMP by TIMP, as well as the prescribed TGF β boundary condition. Compared to one-way coupling, two-way coupling results in elevated measures of stress near the abluminal surface - characteristic of a maladaptive positive feedback loop (Fig. 8( c ) ). These subtle differences may affect the structural remodeling and integrity of the aneurysm. The role of diffusive signaling modeled by PDEs is highlighted in each test problem. For the tissue graft, increased species diffusivity leads to uniform growth from localized morphogen production. In the aneurysm case, diffusivity may be a factor that influences whether a vessel progresses to an aneurysm or simply thickens. This presents an alternate means of restoring homeostasis when compared to two-way coupling, whereby altering the diffusive properties of tissues and chemical species can reduce the heterogeneity of growth and stress distributions. The aneurysm example is interesting because MMP production peaks at a region of half maximal TGF β . When TGF β levels rise above the optimal value, the aneurysm began to split (as seen in case 3). We do not interpret the results of case 3 as a mechanism for the clinical presentation of double aneurysms since the biochemical systems biology is simplified for illustrative purposes. However, our results exemplify how kinematic growth theories may be used to identify potential mechanisms for growth phenomena to be further evaluated experimentally. One theoretical oversight of the aneurysm example is that we simulate growth in response to MMP without remodeling of the matrix constituents (i.e., the material parameters and constituent densities do not change despite implied changes in collagen, elastin, etc.). Nonetheless, the results of aneurysm cases 1 and 2 qualitatively reproduce stress gradients observed in previous works where aneurysm growth was stress-driven without any consideration of the underlying systems biology [ 57 ]. Specifically, we find that the ends of the aneurysm (near axial positions Z = 0.5–1.5) are characterized by high stress gradients with the global maximum value spatially associated toward the healthy vessel and the global minimum value associated with the center of the aneurysm. Additionally, aneurysm case 2 demonstrates slightly greater growth in regions 1 mm from the center of the aneurysm (Fig. 10 ). Similar trends in growth are observed in Ref. 30 ] for cases where vessels experienced axial displacement or cases where the collagen fibrils were oriented > 40 deg diagonally with respect to the axial direction. Previously published studies of G&R have primarily utilized commercial software (e.g., abaqus with user-defined subroutines defined via increasingly outdated FORTRAN code). Additionally, the linearization and discretization of these models are theoretically complicated which belies the inherent simplicity of kinematic growth theory. Collectively, this may hinder reproducibility and adoption of these simulation frameworks for kinematic growth theories since research teams often comprise multidisciplinary biomedical scientists rather than computational physicists with numerical expertise. Through the course of this study, we found the autodifferentiation within fenics to be a straightforward, elegant solution to this hurdle which has motivated us to further pursue the development of autodifferentiation in febio . This would greatly simplify the extension the current G&R computational framework to be dependent on any user-defined variable of interest (e.g., stress, strain, and chemical species activity) and in more physiologically relevant constitutive frameworks (e.g., biphasic poroelasticity) without extensive additional software development. 5 Conclusion We present kinematic growth theory coupled with systems biology in the open-source finite element software febio and fenics . Tissue growth was rate-driven in response to chemical species evolution rates as opposed to tissue stress or strain. Illustrative examples highlight the roles of mechanosensitive coupling and biochemical diffusivity during growth. We also cross-verify two open-source finite element software with nearly identical results. We hope that our freely available open-source software will enable others to apply and expand these developments in the broad field of systems mechanobiology. Data Availability Statement The software and input files for all test problems are freely available and open-source. febio may be downloaded as an executable to develop models in a graphical interface at the following link. 3 The source code and the input files for the febio implementation are available at the following link. The code is included as a solver plugin, which is a dynamic-linked library (.dll) that the user includes at run-time. The source code and input files for the fenics implementation are available at the following link. Supplementary Material Supplementary Material Supplementary PDF Supplementary_Material_BIO-24-1242.zip (8.6MB, zip) Acknowledgment We are thankful to Dr. Gerard Ateshian for discussions related to this work. Appendix A.1 Modifications to the Multiphasic Material Framework in febio . Biological tissues are hydrated materials containing a solid phase and an interstitial fluid phase. Chemical species may move with either the solid (solid-bound molecules) or the fluid phase (solutes). febio simulates realistic biological tissues in a multiphasic framework which can analyze solid, fluid, strain-dependent hydraulic permeability, mechanics, as well as solute and solid-bound molecule transport, reaction kinetics, osmotic pressures, and electrochemical potentials. However, due to the difference in time-scales between the biophysics of daily living and long-term growth and remodeling, we assume that the solutions may be approximated only from solid mechanics, solute kinetics, and solute diffusivity. Similar simplifications are prevalent in the kinematic growth literature [ 49 , 58 – 61 ]. febio utilizes a monolithic solution approach for primary variables (in contrast, to fenics wherein the user could develop an equivalent monolithic solver or separate solvers for mechanics and reaction/transport). Since the simulations in this work utilized only a subset of the capabilities of the multiphasic solver in febio , we developed a separate, simplified version of the multiphasic solver that considers only displacement increments Δ u and concentration increments Δ c α . This was implemented as a solver plugin within febio that solves (hyper-)elastic deformation and chemical species reaction and diffusion kinetics [ 31 ]. To provide equivalence with the physical assumptions implied by the separate solver approach used for the fenics implementation, notable simplifications of the multiphasic solver include: Negligible influence of the solid-volume fraction φ σ (the volume fraction occupied by the nonreactive solid phase, i.e., φ σ = 0 ). This limits the ability of our models to predict solid remodeling where hyperelastic potentials are related to the density of solid constituents (e.g., collagen). Negligible saturation (volume fraction φ α ) of chemical species w.r.t the solvent volume fraction φ w (i.e., ∑ α φ α ≪ φ w ). This eliminates the need to consider viscous drag forces between solutes so that solute transport is analogous to Fick's law of diffusion. It also eliminates the need to explicitly account for solvent supplies (water) or permeability/porosity, which may alternatively be modeled with an effective diffusion coefficient. Chemical concentrations are independent of volumetric deformation (i.e., φ α is constant with respect to Δ u aside from the two-way coupling cases). Deformations due to F g are slow enough that solute mass implicitly increases/decreases with tissue volume. Reaction kinetics are insensitive to stress and strain (i.e., no deviation from idealized kinetics). The signaling kinetics for the two-way coupling cases are only stress sensitive from a phenomenological, systems biology perspective (e.g., hypothesizing that cellular TGF β transcription/production is increased by stress/strain). However, more thermodynamically consistent models could also be used [ 33 , 62 , 63 ]. Solute diffusivity is spatially isotropic and constant [ 49 , 60 ]. These assumptions could be relaxed in the future to allow modeling of anisotropic tissues. Electrochemical and osmotic potentials are negligible and solutes behave ideally. These potentials could be re-added to model swelling in highly hydrated tissues such as cartilage. A.2 Useful Relationships and Derivatives. The fourth-order tensor products of two second-order tensors with an orthonormal Cartesian basis (i.e., e i = { e 1 , e 2 , e 3 } ) are defined as ( A ⊗ B ) ijkl = A i j B k l ( A ⊘ B ) ijkl = A i k B j l ( A ○ \ B ) ijkl = A i l B j k ( A ⊙ B ) ijkl = 1 2 [ A ⊘ B + A ○ \ B ] = 1 2 [ A i k B j l + A i l B j k ] Relevant derivatives are presented here ∂ A ∂ A = I ⊘ I ∂ A T ∂ A = I ○ \ I ∂ A − 1 ∂ A = − A − 1 ⊘ A − T , ∂ A − T ∂ A = − A − T ○ \ A − 1 ∂ tr ( A ) ∂ A = I , ∂ det ( A ) ∂ A = det ( A ) A − T A.3 Governing Equations and Virtual Work. The conservation of momentum for the solid is given by [ 41 ] ρ D v D t = div σ (A1) The conservation of mass for each chemical species α is given as ∂ c α ∂ t + div ( j α ) = ρ ^ α (A2) c α = ρ α M α , c ^ α = ρ ^ α M α (A3) D D t c α + div j α = ∑ i c ^ i α (A4) where M α is the molar mass of species α and i is summed over all reactions. The linearization and discretization of the problems begin with defining the weak form of the problem via the statement of virtual work. The statement of virtual work arising from the axioms of conservation is δ W = δ G + ∫ b δ v · div σ d v + ∑ α ∫ b δ c α [ D D t c α + div j α ] d v (A5) where δ c α is the virtual molar energy of chemical species α . The contribution of reactions to virtual work is given by δ G = ∑ α ∑ i ∫ b δ c α c ^ i α d v (A6) with i is summed over all reactions and α summed over all chemical species. After applying the divergence theorem we separate the statement of virtual work into internal and external contributions given by δ W int = δ G + ∫ b σ : δ d s d v + ∑ α ∫ b [ j α · grad δ c α − δ c α ∂ c α ∂ t ] d v (A7) δ W ext = ∫ ∂ b [ δ vt · n + ∑ α δ c α j α · n ] d a (A8) where δ v is the virtual velocity of the solid phase, d s = ( grad δ v + grad T δ v ) / 2 (not to be confused with diffusivity tensors d α ), b is the body (volume) of interest, n is a normal to the surface of b . A.4 Linearization and Discretization. Although the same statement of virtual work applies to models in febio and fenics , the implementation strategies used to linearize and discretize the solutions differ. In febio , we employed a single solver to a monolithic problem where the linearization and discretization are explicitly derived. In fenics , we employed separate solvers to evaluate nodal displacements and chemical concentrations where the linearization and discretization were achieved through autodifferentiation functions. The linearization and discretization are presented using the notation corresponding to each software. A.4.1 Linearization in febio . The incremental-iterative solution strategy in febio emanates from linearization of the virtual work about a known configuration. Given the assumption of independence between deformation and concentration incorporated in this work, linearization of the above virtual work statement can be written symbolically as δ W + D δ W [ Δ u ] + ∑ α D δ W [ Δ c α ] ≈ 0 (A9) where D δ W [ • ] indicates the directional derivative of δ W in the direction of • ( Δ u or Δ c α ). Linearization along a displacement increment Δ u is given by D ( S : δ E ˙ ) [ Δ u ] d V = [ δ d s : C : Δ ε + σ : ( grad T Δ u · grad δ v ) ] d v (A10) where E ˙ = F T · D · F is the material time derivative of the strain rate deformation tensor, Δ ε = ( grad Δ u + grad T Δ u ) / 2 is the infinitesimal strain tensor, and C is the fourth-order spatial elasticity tensor given by C = J − 1 [ F ⊘ F ] : C : [ F T ⊘ F T ] (A11) The material elasticity tensor C is composed of elastic contributions from C e and growth contributions C g [ 35 ] C = 2 d S ( F , F g ) d C = C e + C g (A12) The elastic contribution is given by C e = 2 ∂ S ( F , F g ) ∂ C | F g = J g [ F g − 1 ⊘ F g − 1 ] : L e : [ F g − T ⊘ F g − T ] (A13) L e = 2 ∂ S e ∂ C e | F g (A14) The growth contribution is given from C g = 2 [ ∂ S ( F , F g ) ∂ F g | F : ∂ F g ∂ ϑ ] ⊗ ∂ ϑ ∂ C (A15) The first term is given by [ 35 ] ∂ S ∂ F g = S ⊗ F g − T − [ F g − 1 ⊘ S + S ○ \ F g − 1 ] − J g [ F g − 1 ⊘ F g − 1 ] : 1 2 L e : [ F g − T ○ \ C e + C e ⊘ F g − T ] (A16) The second term is determined by the type of growth. Symmetric linear and area type growth have the respective terms: ∂ F g ∂ ϑ = { n ∥ ⊗ n ∥ line − type growth 1 2 ϑ [ I − n ⊥ ⊗ n ⊥ ] area − type growth (A17) where n ∥ is along the direction of growth and n ⊥ is in the direction perpendicular to the plane of area growth. The final term is found by evaluating the residual R and its tangent K : R = ϑ − ϑ n − ϑ ˙ Δ t (A18) K = ∂ R ∂ ϑ = 1 − Δ t ∂ ϑ ˙ ∂ ϑ (A19) ∂ ϑ ˙ ∂ ϑ = ϕ ( c ˙ ) ∂ k ( ϑ ) ∂ ϑ + k ( ϑ ) ∂ ϕ ( c ˙ ) ∂ ϑ (A20) The function ϕ ( c ˙ ) is independent of ϑ ; thus ∂ ϑ ˙ ∂ ϑ = ϕ ( c ˙ ) ∂ k ( ϑ ) ∂ ϑ (A21) K = 1 − Δ t ϕ ( c ˙ ) ∂ k ( ϑ ) ∂ ϑ (A22) The residual is evaluated with respect to C and rearranged to solve for ∂ ϑ ∂ C ∂ R ∂ C = ∂ ϑ ∂ C − ∂ ϑ ˙ ∂ C Δ t = 0 (A23) ∂ ϑ ∂ C = ∂ ϑ ˙ ∂ C Δ t = [ ϕ ( c ˙ α ) ∂ k ( ϑ ) ∂ ϑ ∂ ϑ ∂ C + k ( ϑ ) ∂ ϕ ( c ˙ α ) ∂ c ˙ α ∂ c ˙ α ∂ C ] Δ t (A24) ∂ ϑ ∂ C = k ( ϑ ) Δ t K ∂ ϕ ( c ˙ α ) ∂ c ˙ α ∂ c ˙ α ∂ C (A25) where ∂ ϕ ( c ˙ α ) ∂ c ˙ α = { ∂ ϕ ( c ˙ α ) ∂ c ˙ α = growth species 0 α ≠ growth species (A26) ∂ k ( ϑ ) ∂ ϑ = ϑ max γ ( exp ( − ϑ + ϑ a γ ) ( 1 + exp ( − ϑ + ϑ a γ ) ) 2 + exp ( − ϑ − ϑ a γ ) ( 1 + exp ( − ϑ − ϑ a γ ) ) 2 ) (A27) ∂ c ˙ α ∂ C = J 2 F − 1 · c ^ ε α · F − T (A28) c ^ ε α = ∂ c ^ α ∂ ε = ∑ i ν i α ζ ^ ε , i α (A29) ζ ^ ε , i = ∂ ζ ^ i ∂ ε = { 0 one − way coupling ∂ ζ ^ i ∂ tr ( σ ) I : C two − way coupling (A30) The next two linearization terms in the direction of Δ u are 0 due to the assumptions made about minimal interactions between mechanics and chemical kinetics due to disparate timescales, i.e., D ( J α · Grad δ c α ) [ Δ u ] d V = 0 , D ( δ c α ∂ c α ∂ t ) [ Δ u ] d V = 0 (A31) where J α = J F − 1 j α F − T . The final linearization term in the direction of Δ u is given by D ( δ G ) [ Δ u ] d V = ∑ α ∫ b δ c α [ c ^ ε α : Δ ε ] d v = ∑ α ∑ i ν i α ∫ b δ c α [ ζ ^ ε , i : Δ ε ] d v (A32) Linearization along concentration increment Δ c γ are given by D ( S : δ E ˙ ) [ Δ c γ ] d V = Δ c γ ( σ γ ′ : δ d s ) d v , (A33) D ( J α · Grad δ c α ) [ Δ c γ ] d V = δ α γ grad δ c γ · d α · grad Δ c γ d v (A34) D ( ∂ c α ∂ t δ c α ) [ Δ c γ ] d V = δ c α δ α γ Δ t Δ c γ d v (A35) D δ G [ Δ c γ ] d V = ∑ α ∑ i ∫ b δ c α ∂ c ^ i α ∂ c γ Δ c γ d v (A36) Here, we include the functions σ γ ′ = J − 1 F · ∂ S ∂ c γ · F T (A37) ∂ S ∂ c γ = ∂ S ∂ F g : ∂ F g ∂ ϑ ∂ ϑ ∂ c γ (A38) ∂ ϑ ∂ c γ = ∂ ϑ ∂ ϑ ˙ ∂ ϑ ˙ ∂ c ˙ γ ∂ c ˙ γ ∂ c γ = k ( ϑ ) ∂ ϕ ( c ˙ ) ∂ c γ (A39) Lastly, we linearize chemical reactions for Eq. (A36) . Mass-action forward molar production rates for a species α are linearized with respect to a species β according to ∂ c ^ i α ∂ c β = ζ ^ i ( ν i α ) 2 c α (A40) For Hill-type reaction kinetics (Eq. (14) ), the molar supply (Eq. (13) ) is linearized as ∂ c ^ → β α ∂ c β = k H ∂ f → β α ∂ c β . The linearization for Hill functions f → β α where a species α is activated in response to a species β are given as ∂ f → β α ∂ c β = n K n f → β α c β ( K n + ( c β ) n ) (A41) The tangent of the activation reaction that occurs due to species α is thus ∂ c ^ → β α ∂ c β = k H ∂ f → β α ∂ c β (A42) The tangent of the inhibition reaction f ¬ γ α that occurs due to species γ is similarly given as ∂ c ^ ¬ γ α ∂ c γ = ∂ ( k H ( 1 − f → γ α ) ) ∂ c γ = − k H ∂ f → γ α ∂ c γ (A43) The tangents of regulation with activation of α by β and inhibition by γ given by with respect to concentration are given by ∂ c ^ → β ∧ ¬ γ α ∂ c β = k H f ¬ γ α ∂ f → β α ∂ c α (A44) ∂ c ^ → β ∧ ¬ γ α ∂ c γ = − k H f → β α ∂ f → γ α ∂ c β (A45) A.4.2 Derivation of Linearization for External Work in febio . The linearization of external virtual work is given by D ( δ W ext e ) [ Δ u ] = ∫ ∂ b ( t · n δ v + ∑ α δ c α j α · n ) · ( ∂ Δ u ∂ η 1 × g 2 + g 1 × ∂ Δ u ∂ η 2 ) d η 1 d η 2 (A46) with the elemental area on ∂ b given by d a = | g 1 × g 2 | d η 1 d η 2 where η 1 , η 2 are the parametric coordinates of the element, g i = ∂ x / ∂ η i , and the outward normal vector is n = ( g 1 × g 2 ) / | g 1 × g 2 | . A.4.3 Derivation of Discretization in febio . Virtual and nodal variables are interpolated through the shape functions N a and N b as δ v = ∑ a N a δ v a , Δ u = ∑ b N b Δ , u b (A47) δ c α = ∑ a N a δ c a α , Δ c γ = ∑ b N b Δ c b γ (A48) We discretize δ W int e as δ W int e = ∑ e = 1 n e ∑ k = 1 n int e W k J η ∑ a = 1 m [ δ v δ c a α … δ c a γ ] · [ r a u r a α ⋮ r a γ ] (A49) r a u = σ · grad N a (A50) r a α = j α · grad N a − N a ∂ c α ∂ t + N a ∑ i c ^ i α (A51) We discretize D δ W int = D δ W int [ Δ u ] + ∑ γ D δ W int [ Δ c γ ] as D δ W int = ∑ e = 1 n e ∑ k = 1 n int e W k J η ∑ a = 1 m ∑ b = 1 m [ δ v a δ c a α ⋯ δ c a γ ] · [ K a b u u k a b u α ⋯ k a b u γ k a b α u k a b α α ⋯ k a b α γ ⋮ ⋮ ⋱ ⋮ k a b γ u k a b γ α ⋯ k a b γ γ ] · [ Δ u Δ c α ⋮ Δ c γ ] (A52) where n e is the number of elements, n int e is the number of elemental integration points, m is the number of nodes in the element, W k is the Gauss weight of the k th integration point, and J η = | g 1 × g 2 | . The internal stiffness matrix components are given as K a b u u = grad N a · C e · grad N b + ( grad N a · σ · grad N b ) I (A53) k a b α u = N a ∑ i c ^ ε , i α · grad N b (A54) k a b u α = N b σ α ′ · grad N a (A55) k a b α γ = − δ α γ grad N a · d α · grad N b − N a N b 1 Δ t + N a N b ∑ α ∑ i ∂ c ^ i α ∂ c γ (A56) We discretize δ W ext − D δ W ext = 0 as δ W ext = ∑ e = 1 n e ∑ k = 1 n int e W k J η ∑ a = 1 m ∑ b = 1 m [ δ v a δ c a α ⋯ δ c a γ ] · [ N a t N a ( j α · n ) ⋮ N a ( j γ · n ) ] (A57) − D δ W ext = ∑ e = 1 n e ∑ k = 1 n int e W k J η ∑ a = 1 m ∑ b = 1 m [ δ v a δ c a α ⋯ δ c a γ ] · [ K a b u u 0 ⋯ 0 k a b α u 0 ⋯ 0 ⋮ ⋮ ⋱ ⋮ k a b γ u 0 ⋯ 0 ] · [ Δ u Δ c α ⋮ Δ c γ ] (A58) Here, we define the external stiffness matrix components K a b u u = ( t · n ) N a A { ∂ N b ∂ η 1 g 2 − ∂ N a ∂ η 2 g 1 } (A59) k a b α u = − ( j α · n ) N a { ∂ N b ∂ η 1 g 2 − ∂ N a ∂ η 2 g 1 } × n (A60) where A is the antisymmetric tensor with dual vector v such that A { v } = v × q . A.4.4 Linearization and Discretization in fenics . We opted to enforce conservation of linear momentum and conservation of mass using separate solvers in fenics because of the disparate time-scales for mechanics and chemical kinetics in the problems we considered. Because of this assumption, linearization of the separate statements of virtual work is greatly simplified. In fenics , the linearization of each virtual work statement is calculated via the built-in autodifferentiation functions supplied by the fenics form language (FFL). The user specifies each statement of virtual work along with the functional dependence on primary variables. The variational form (statement of virtual work) for hyperelastic displacements was stated in relation to the first Piola–Kirchhoff stress P = J σ F − T F ( u , δ u ) ≡ δ W e ( u , δ u ) = ∫ ( P : Grad δ v ) d V + ( J F − T · t ) · δ v d a (A61) where ( J F − T · t ) is the pull-back of external surface pressures (note: this is energetically conjugate with F in a different intermediate configuration than the “fictitious” intermediate configuration posited by kinematic growth theory). The Gâteaux derivative of F ( u , δ u ) is given by J F ( u , δ u ) = D δ W e ( u , δ u ) [ Δ u ] = d d ϵ δ W e ( u + ϵ Δ u | ϵ = 0 ) (A62) which is solved by fenics 's built-in autodifferentiation function. The nonlinear solver is assembled by providing F ( u , δ u ) , J F ( u , δ u ) , the boundary conditions, and an initial guess of the solution to fenics . The variational form for the molar energy (conservation of solute mass) is given by F ( c , δ c ) ≡ δ W e ( c , δ c ) = ∫ b ∑ α ( Δ t − 1 ( c n + 1 α − c n α ) d v + ∫ b ∑ α ( D α ∇ c n + 1 α · ∇ δ c n α ) d v − ∫ b ∑ i ( c ^ i δ c n α ) d v (A63) The Gâteaux derivative of F ( c , δ c ) is similarly given by J F ( c , δ c ) = D δ W e ( c , δ c ) [ Δ c ] = d d ϵ δ W e ( c + ϵ Δ c | ϵ = 0 ) (A64) which is also solved by autodifferentiation. A separate nonlinear solver is assembled by providing F ( c , δ c ) , J F ( c , δ c ) , the boundary conditions, and an initial guess of the solution to fenics . Footnotes 3 https://febio.org 4 https://fenicsproject.org 5 https://github.com/febiosoftware/FEBioERD 6 https://github.com/amir-cardiolab/Multiscale_mechanobiology_GR Funding Data NSF OAC (Award No. 2246911; Funder ID: 10.13039/100000001). NIH (Award Nos. 2R01GM083925, 1R01CA290182, and 5R01CA290182; Funder ID: 10.13039/100000002). References [1]. Ambrosi, D., Ateshian, G. A., Arruda, E. M., Cowin, S. C., Dumais, J., Goriely, A., Holzapfel, G. A., et al., 2011, “Perspectives on Biological Growth and Remodeling,” J. Mech. Phys. Solids, 59(4), pp. 863–883. 10.1016/j.jmps.2010.12.011 [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [2]. Irons, L., Latorre, M., and Humphrey, J. D., 2021, “From Transcript to Tissue: Multiscale Modeling From Cell Signaling to Matrix Remodeling,” Ann. Biomed. Eng., 49(7), pp. 1701–1715. 10.1007/s10439-020-02713-8 [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [3]. Carlier, A., Geris, L., Lammens, J., and Van Oosterwyck, H., 2015, “Bringing Computational Models of Bone Regeneration to the Clinic,” WIREs Syst. Biol. Med., 7(4), pp. 183–194. 10.1002/wsbm.1299 [ DOI ] [ PubMed ] [ Google Scholar ] [4]. Carlier, A., Lammens, J., Van Oosterwyck, H., and Geris, L., 2015, “Computational Modeling of Bone Fracture Non-Unions: Four Clinically Relevant Case Studies,” Silico Cell Tissue Sci., 2(1), pp. 1–12. 10.1186/s40482-015-0004-x [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [5]. Guo, Y., Mofrad, M. R. K., and Tepole, A. B., 2022, “On Modeling the Multiscale Mechanobiology of Soft Tissues: Challenges and Progress,” Biophys. Rev., 3(3), p. 031303. 10.1063/5.0085025 [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [6]. Skalak, R., Dasgupta, G., Moss, M., Otten, E., Dullemeijer, P., and Vilmann, H., 1982, “Analytical Description of Growth,” J. Theor. Biol., 94(3), pp. 555–577. 10.1016/0022-5193(82)90301-0 [ DOI ] [ PubMed ] [ Google Scholar ] [7]. Skalak, R., 1981, Growth as a Finite Displacement Field, Springer, The Netherlands, pp. 347–355. [ Google Scholar ] [8]. Rodriguez, E. K., Hoger, A., and McCulloch, A. D., 1994, “Stress-Dependent Finite Growth in Soft Elastic Tissues,” J. Biomech., 27(4), pp. 455–467. 10.1016/0021-9290(94)90021-3 [ DOI ] [ PubMed ] [ Google Scholar ] [9]. Humphrey, J., and Rajagopal, K., 2002, “A Constrained Mixture Model for Growth and Remodeling of Soft Tissues,” Math. Models Methods Appl. Sci., 12(3), pp. 407–430. 10.1142/S0218202502001714 [ DOI ] [ Google Scholar ] [10]. Ateshian, G. A., and Humphrey, J. D., 2012, “Continuum Mixture Models of Biological Growth and Remodeling: Past Successes and Future Opportunities,” Annu. Rev. Biomed. Eng., 14(1), pp. 97–111. 10.1146/annurev-bioeng-071910-124726 [ DOI ] [ PubMed ] [ Google Scholar ] [11]. Ashrafi, M., Gubaua, J. E., Pereira, J. T., Gahlichi, F., and Doblaré, M., 2020, “A Mechano-Chemo-Biological Model for Bone Remodeling With a New Mechano-Chemo-Transduction Approach,” Biomech. Model. Mechanobiol., 19(6), pp. 2499–2523. 10.1007/s10237-020-01353-0 [ DOI ] [ PubMed ] [ Google Scholar ] [12]. Irons, L., and Humphrey, J. D., 2020, “Cell Signaling Model for Arterial Mechanobiology,” PLoS Comput. Biol., 16(8), p. e1008161. 10.1371/journal.pcbi.1008161 [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [13]. Kim, J., and Wagenseil, J. E., 2015, “Bio-Chemo-Mechanical Models of Vascular Mechanics,” Ann. Biomed. Eng., 43(7), pp. 1477–1487. 10.1007/s10439-014-1201-7 [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [14]. Sadrabadi, M. S., Eskandari, M., Feigenbaum, H. P., and Arzani, A., 2021, “Local and Global Growth and Remodeling in Calcific Aortic Valve Disease and Aging,” J. Biomech., 128,p. 110773. 10.1016/j.jbiomech.2021.110773 [ DOI ] [ PubMed ] [ Google Scholar ] [15]. Zeigler, A., Richardson, W., Holmes, J., and Saucerman, J., 2016, “A Computational Model of Cardiac Fibroblast Signaling Predicts Context-Dependent Drivers of Myofibroblast Differentiation,” J. Mol. Cell. Cardiol., 94, pp. 72–81. 10.1016/j.yjmcc.2016.03.008 [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [16]. Humphrey, J. D., Dufresne, E. R., and Schwartz, M. A., 2014, “Mechanotransduction and Extracellular Matrix Homeostasis,” Nat. Rev. Mol. Cell Biol., 15(12), pp. 802–812. 10.1038/nrm3896 [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [17]. Loerakker, S., and Ristori, T., 2020, “Computational Modeling for Cardiovascular Tissue Engineering: The Importance of Including Cell Behavior in Growth and Remodeling Algorithms,” Curr. Opin. Biomed. Eng., 15, pp. 1–9. 10.1016/j.cobme.2019.12.007 [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [18]. Sree, V. D., and Tepole, A. B., 2020, “Computational Systems Mechanobiology of Growth and Remodeling: Integration of Tissue Mechanics and Cell Regulatory Network Dynamics,” Curr. Opin. Biomed. Eng., 15, pp. 75–80. 10.1016/j.cobme.2020.01.002 [ DOI ] [ Google Scholar ] [19]. Hahn, C., and Schwartz, M. A., 2009, “Mechanotransduction in Vascular Physiology and Atherogenesis,” Nat. Rev. Mol. Cell Biol., 10(1), pp. 53–62. 10.1038/nrm2596 [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [20]. Freedman, B. R., Bade, N. D., Riggin, C. N., Zhang, S., Haines, P. G., Ong, K. L., and Janmey, P. A., 2015, “The (Dys) Functional Extracellular Matrix,” Biochim. Biophys. Acta (BBA)-Mol. Cell Res., 1853(11 Pt B), pp. 3153–3164. 10.1016/j.bbamcr.2015.04.015 [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [21]. Sagi, I., and Gaffney, J., 2015, Matrix Metalloproteinase Biology, Wiley, Hoboken, NJ. [ Google Scholar ] [22]. Sugita, S., and Matsumoto, T., 2017, “Multiphoton Microscopy Observations of 3D Elastin and Collagen Fiber Microstructure Changes During Pressurization in Aortic Media,” Biomech. Model. Mechanobiol., 16(3), pp. 763–773. 10.1007/s10237-016-0851-9 [ DOI ] [ PubMed ] [ Google Scholar ] [23]. Yoshida, K., McCulloch, A. D., Omens, J. H., and Holmes, J. W., 2020, “Predictions of Hypertrophy and Its Regression in Response to Pressure Overload,” Biomech. Model. Mechanobiol., 19(3), pp. 1079–1089. 10.1007/s10237-019-01271-w [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [24]. Hayenga, H. N., Thorne, B. C., Peirce, S. M., and Humphrey, J. D., 2011, “Ensuring Congruency in Multiscale Modeling: Towards Linking Agent Based and Continuum Biomechanical Models of Arterial Adaptation,” Ann. Biomed. Eng., 39(11), pp. 2669–2682. 10.1007/s10439-011-0363-9 [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [25]. Pivonka, P., Zimak, J., Smith, D. W., Gardiner, B. S., Dunstan, C. R., Sims, N. A., Martin, T. J., and Mundy, G. R., 2008, “Model Structure and Control of Bone Remodeling: A Theoretical Study,” Bone, 43(2), pp. 249–263. 10.1016/j.bone.2008.03.025 [ DOI ] [ PubMed ] [ Google Scholar ] [26]. Rouhi, G., Epstein, M., Sudak, L., and Herzog, W., 2007, “Modeling Bone Resorption Using Mixture Theory With Chemical Reactions,” J. Mech. Mater. Struct., 2(6), pp. 1141–1155. 10.2140/jomms.2007.2.1141 [ DOI ] [ Google Scholar ] [27]. Schäfer, A., Weickenmeier, J., and Kuhl, E., 2019, “The Interplay of Biochemical and Biomechanical Degeneration in Alzheimer's Disease,” Comput. Methods Appl. Mech. Eng., 352, pp. 369–388. 10.1016/j.cma.2019.04.028 [ DOI ] [ Google Scholar ] [28]. Peirce-Cottler, S. M., Sander, E. A., Fisher, M. B., Deymier, A. C., LaDisa, J. F. Jr., O'Connell, G., Corr, D. T., et al., 2024, “A Systems Approach to Biomechanics, Mechanobiology, and Biotransport,” ASME J. Biomech. Eng., 146(4), p. 040801. 10.1115/1.4064547 [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [29]. Maas, S., Ellis, B. J., Ateshian, G. A., and Weiss, J. A., 2012, “Febio: Finite Elements for Biomechanics,” ASME J. Biomech. Eng., 134(1), p. 011005. 10.1115/1.4005694 [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [30]. Logg, A., Wells, G. N., and Hake, J., 2012, “DOLFIN: A C++/Python Finite Element Library,” Computational Science and Engineering (Lecture Notes in chapter 10), Vol. 84, Springer, Berlin Heidelberg, pp. 173–225. [ Google Scholar ] [31]. Maas, S. A., LaBelle, S. A., Ateshian, G. A., and Weiss, J. A., 2018, “A Plugin Framework for Extending the Simulation Capabilities of Febio,” Biophys. J., 115(9), pp. 1630–1637. 10.1016/j.bpj.2018.09.016 [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [32]. Estrada, A. C., Yoshida, K., Saucerman, J. J., and Holmes, J. W., 2021, “A Multiscale Model of Cardiac Concentric Hypertrophy Incorporating Both Mechanical and Hormonal Drivers of Growth,” Biomech. Model. Mechanobiol., 20(1), pp. 293–307. 10.1007/s10237-020-01385-6 [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [33]. Senthilnathan, C., 2024, “Understanding the Mechanics of Growth: A Large Deformation Theory for Coupled Swelling-Growth and Morphogenesis of Soft Biological Systems,” Ph.D. thesis, Massachusetts Institute of Technology, Cambridge, MA. [ Google Scholar ] [34]. Eskandari, M., and Kuhl, E., 2015, “Systems Biology and Mechanics of Growth,” Wiley Interdiscip. Rev. Syst. Biol. Med., 7(6), pp. 401–412. 10.1002/wsbm.1312 [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [35]. Göktepe, S., Abilez, O. J., Parker, K. K., and Kuhl, E., 2010, “A Multiscale Model for Eccentric and Concentric Cardiac Growth Through Sarcomerogenesis,” J. Theor. Biol., 265(3), pp. 433–442. 10.1016/j.jtbi.2010.04.023 [ DOI ] [ PubMed ] [ Google Scholar ] [36]. Menzel, A., 2005, “Modelling of Anisotropic Growth in Biological Tissues,” Biomech. Model. Mechanobiol., 3(3), pp. 147–171. 10.1007/s10237-004-0047-6 [ DOI ] [ PubMed ] [ Google Scholar ] [37]. Göktepe, S., Abilez, O. J., and Kuhl, E., 2010, “A Generic Approach Towards Finite Growth With Examples of Athlete's Heart, Cardiac Dilation, and Cardiac Wall Thickening,” J. Mech. Phys. Solids, 58(10), pp. 1661–1680. [ Google Scholar ] [38]. Rausch, M., Dam, A., Göktepe, S., Abilez, O., and Kuhl, E., 2011, “Computational Modeling of Growth: Systemic and Pulmonary Hypertension in the Heart,” Biomech. Model. Mechanobiol., 10(6), pp. 799–811. 10.1007/s10237-010-0275-x [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [39]. Ambrosi, D., Ben Amar, M., Cyron, C. J., DeSimone, A., Goriely, A., Humphrey, J. D., and Kuhl, E., 2019, “Growth and Remodelling of Living Tissues: Perspectives, Challenges and Opportunities,” J. R. Soc. Interface, 16(157), p. 20190233. 10.1098/rsif.2019.0233 [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [40]. Ateshian, G. A., LaBelle, S. A., and Weiss, J. A., 2024, “Continuum Growth Mechanics: Reconciling Two Common Frameworks,” ASME J. Biomech. Eng., 146(10), p. 101003. 10.1115/1.4065309 [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [41]. Ateshian, G. A., Nims, R. J., Maas, S. A., and Weiss, J. A., 2014, “Computational Modeling of Chemical Reactions and Interstitial Growth and Remodeling Involving Charged Solutes and Solid-Bound Molecules,” Biomech. Model. Mechanobiol., 13(5), pp. 1105–1120. 10.1007/s10237-014-0560-1 [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [42]. Prud'homme, R., 2010, “Flows Reactive Fluids,” 94, Springer, New York. [ Google Scholar ] [43]. Menzel, A., and Kuhl, E., 2012, “Frontiers in Growth and Remodeling,” Mech. Res. Commun., 42, pp. 1–14. 10.1016/j.mechrescom.2012.02.007 [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [44]. Cyron, C. J., Aydin, R. C., and Humphrey, J. D., 2016, “A Homogenized Constrained Mixture (and Mechanical Analog) Model for Growth and Remodeling of Soft Tissue,” Biomech. Model. Mechanobiol., 15(6), pp. 1389–1403. 10.1007/s10237-016-0770-9 [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [45]. Kuhl, E., Maas, R., Himpel, G., and Menzel, A., 2007, “Computational Modeling of Arterial Wall Growth. Attempts Towards Patient-Specific Simulations Based on Computer Tomography,” Biomech. Model. Mechanobiol., 6(5), pp. 321–331. 10.1007/s10237-006-0062-x [ DOI ] [ PubMed ] [ Google Scholar ] [46]. Sheth, R. A., Maricevich, M., and Mahmood, U., 2010, “In Vivo Optical Molecular Imaging of Matrix Metalloproteinase Activity in Abdominal Aortic Aneurysms Correlates With Treatment Effects on Growth Rate,” Atherosclerosis, 212(1), pp. 181–187. 10.1016/j.atherosclerosis.2010.05.012 [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [47]. Gierig, M., Wriggers, P., and Marino, M., 2021, “Computational Model of Damage-Induced Growth in Soft Biological Tissues Considering the Mechanobiology of Healing,” Biomech. Model. Mechanobiol., 20(4), pp. 1297–1315. 10.1007/s10237-021-01445-5 [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [48]. Marino, M., Pontrelli, G., Vairo, G., and Wriggers, P., 2017, “Chapter 4 - Coupling Microscale Transport and Tissue Mechanics: Modeling Strategies for Arterial Multiphysics,” Modeling of Microscale Transport in Biological Processes, Becker S. M., ed., Academic Press, Cambridge, MA, pp. 77–112. [ Google Scholar ] [49]. Marino, M., Pontrelli, G., Vairo, G., and Wriggers, P., 2017, “A Chemo-Mechano-Biological Formulation for the Effects of Biochemical Alterations on Arterial Mechanics: The Role of Molecular Transport and Multiscale Tissue Remodelling,” J. R. Soc. Interface, 14(136), p. 20170615. 10.1098/rsif.2017.0615 [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [50]. Tan, P. M., Buchholz, K. S., Omens, J. H., McCulloch, A. D., and Saucerman, J. J., 2017, “Predictive Model Identifies Key Network Regulators of Cardiomyocyte Mechano-Signaling,” PLoS Comput. Biol., 13(11), p. e1005854. 10.1371/journal.pcbi.1005854 [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [51]. Taber, L. A., 1995, “Biomechanics of Growth, Remodeling, and Morphogenesis,” ASME Appl. Mech. Rev., 48(8), pp. 487–545. 10.1115/1.3005109 [ DOI ] [ Google Scholar ] [52]. Kerckhoffs, R. C. P., Omens, J. H., and McCulloch, A. D., 2012, “A Single Strain-Based Growth Law Predicts Concentric and Eccentric Cardiac Growth During Pressure and Volume Overload,” Mech. Res. Commun., 42, pp. 40–50. 10.1016/j.mechrescom.2011.11.004 [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [53]. Schwarz, E. L., Pfaller, M. R., Szafron, J. M., Latorre, M., Lindsey, S. E., Breuer, C. K., Humphrey, J. D., and Marsden, A. L., 2023, “A Fluid–Solid-Growth Solver for Cardiovascular Modeling,” Comput. Methods Appl. Mech. Eng., 417, p. 116312. 10.1016/j.cma.2023.116312 [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [54]. Teixeira, F. S., Neufeld, E., Kuster, N., and Watton, P. N., 2020, “Modeling Intracranial Aneurysm Stability and Growth: An Integrative Mechanobiological Framework for Clinical Cases,” Biomech. Model. Mechanobiol., 19(6), pp. 2413–2431. 10.1007/s10237-020-01351-2 [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [55]. Keshavarzian, M., Meyer, C. A., and Hayenga, H. N., 2018, “Mechanobiological Model of Arterial Growth and Remodeling,” Biomech. Model. Mechanobiol., 17(1), pp. 87–101. 10.1007/s10237-017-0946-y [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [56]. Scheiner, S., Pivonka, P., and Hellmich, C., 2013, “Coupling Systems Biology With Multiscale Mechanics, for Computer Simulations of Bone Remodeling,” Comput. Methods Appl. Mech. Eng., 254, pp. 181–196. 10.1016/j.cma.2012.10.015 [ DOI ] [ Google Scholar ] [57]. Latorre, M., and Humphrey, J. D., 2020, “Fast, Rate-Independent, Finite Element Implementation of a 3D Constrained Mixture Model of Soft Tissue Growth and Remodeling,” Comput. Methods Appl. Mech. Eng., 368, p. 113156. 10.1016/j.cma.2020.113156 [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [58]. Ambrosi, D., Pezzuto, S., Riccobelli, D., Stylianopoulos, T., and Ciarletta, P., 2017, “Solid Tumors Are Poroelastic Solids With a Chemo-Mechanical Feedback on Growth,” J. Elast., 129(1–2), pp. 107–124. 10.1007/s10659-016-9619-9 [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [59]. Paukner, D., Humphrey, J. D., and Cyron, C. J., 2024, “Multiscale Homogenized Constrained Mixture Model of the Bio-Chemo-Mechanics of Soft Tissue Growth and Remodeling,” Biomech. Model. Mechanobiol., 23(6), pp. 2115–2136. 10.1007/s10237-024-01884-w [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] [60]. Soleimani, M., Muthyala, N., Marino, M., and Wriggers, P., 2020, “A Novel Stress-Induced Anisotropic Growth Model Driven by Nutrient Diffusion: Theory, FEM Implementation and Applications in Bio-Mechanical Problems,” J. Mech. Phys. Solids, 144, p. 104097. 10.1016/j.jmps.2020.104097 [ DOI ] [ Google Scholar ] [61]. Verner, S. N., and Garikipati, K., 2018, “A Computational Study of the Mechanisms of Growth-Driven Folding Patterns on Shells, With Application to the Developing Brain,” Ext. Mech. Lett., 18, pp. 58–69. 10.1016/j.eml.2017.11.003 [ DOI ] [ Google Scholar ] [62]. Afshar, A., and Di Leo, C. V., 2021, “A Thermodynamically Consistent Gradient Theory for Diffusion–Reaction–Deformation in Solids: Application to Conversion-Type Electrodes,” J. Mech. Phys. Solids, 151, p. 104368. 10.1016/j.jmps.2021.104368 [ DOI ] [ Google Scholar ] [63]. Chockalingam, S., and Cohen, T., 2024, “A Large Deformation Theory for Coupled Swelling and Growth With Application to Growing Tumors and Bacterial Biofilms,” J. Mech. Phys. Solids, 187, p. 105627. 10.1016/j.jmps.2024.105627 [ DOI ] [ Google Scholar ] Associated Data This section collects any data citations, data availability statements, or supplementary materials included in this article. Supplementary Materials Supplementary Material Supplementary PDF Supplementary_Material_BIO-24-1242.zip (8.6MB, zip) Data Availability Statement The software and input files for all test problems are freely available and open-source. febio may be downloaded as an executable to develop models in a graphical interface at the following link. 3 The source code and the input files for the febio implementation are available at the following link. The code is included as a solver plugin, which is a dynamic-linked library (.dll) that the user includes at run-time. The source code and input files for the fenics implementation are available at the following link. Articles from Journal of Biomechanical Engineering are provided here courtesy of American Society of Mechanical Engineers ACTIONS View on publisher site PDF (3.8 MB) Cite Collections Permalink RESOURCES Similar articles Cited by other articles Links to NCBI Databases Follow NCBI NCBI on X (formerly known as Twitter) NCBI on Facebook NCBI on LinkedIn NCBI on GitHub NCBI RSS feed Connect with NLM NLM on X (formerly known as Twitter) NLM on Facebook NLM on YouTube National Library of Medicine 8600 Rockville Pike Bethesda, MD 20894 Web Policies FOIA HHS Vulnerability Disclosure Help Accessibility Careers NLM NIH HHS USA.gov Back to Top

Related documents

• Mechano-molecular Interplay: A Multiscale Approach to Mechanically Driven Diseases#191812
• Systems and methods for providing modular architectures for robotic end … — Boston Dynamics, Inc. (US20250033225A1)#607444
• Kinematic consistency relations of large-scale structures#1053395
• Geometric mechanogenomics: engineering boundary conditions for deterministic cell fate control.#653458
• Computational Mechanics#253469