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

Abstract

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

Fig. 1

a Time-dependent dissociation of NaCl into Na⁺ and Cl⁻. The ion concentration c(t) was evaluated from a finite element analysis and from the analytical solution of (53), with c_t = 1, k_F = 1, and K_a = 10. b Time discretization error analysis

Fig. 2

a Comparison of finite element results (thick gray curves) against analytical solution (thin black curves) for the cartilage matrix biosynthesis, binding and degradation study by DiMicco and Sah (2003). The function f (α, x/h) represents the normalized spatial distribution of the soluble matrix product concentration, where c^u = k_f f/k_b, and bound matrix concentration, where c^b = k_f f/k_d, over the range 0 ≤ x ≤ h, for different values of the non-dimensional parameter α (see text). b Mesh convergence results for α = 8

Fig. 3

Spatial distribution of the solid apparent density at steady state, in a prismatic bar of length h = 10, subjected to an axial traction σ₀ = −10 and body force per volume f₀ = −1, over the range 0 ≤ x ≤ h, with c = 10⁴, γ = 2, and ψ₀ = 0.01. The finite element solution is obtained using B = 1 in (44), with initial condition ${\rho }_r^s=1$ at t = 0, and with the analysis performed until t = 200. a Prismatic bar geometry and spatial distribution of ${\rho }_r^s$ at steady state. b Comparison of finite element analysis and analytical solution of (55) over the entire range of x. c Mesh convergence results

Fig. 4

a Surface displacement u (h, t)/hξ for 1D analysis of interstitial solid growth. The analytical solution is based on (59), with $\eta =(1-{\phi }_r^s)\widehat{\zeta }\bar{\mathcal{V}}\tau$. The finite element solution was obtained using h = 1, H_A = 1, k = 1, and $(1-{\phi }_r^s)\widehat{\zeta }\bar{\mathcal{V}}={10}^{-3}$. b Mesh convergence results

Fig. 5

Axisymmetric analysis of engineered cartilage growth with glucose consumption and synthesis of negatively charged CS. The tissue construct, which is initially cylindrical, rests on the impermeable bottom of the culture dish. This cross-sectional view shows an outline of the construct and contour maps of the glucose concentration, $c_r^F$, p and ${\phi }_r^s$ at day 30

Fig. 6

Cartilage construct engineered from chondrocyte-seeded agarose gels. Side view of construct at right shows an overlay of a day 0 construct over the image of the day 53 construct. The increased bulging and construct opacity at the outer periphery is qualitatively consistent with the shape and matrix deposition distribution of the finite element results in Fig. 5 (Scale bar 1mm)