A note on stress-driven anisotropic diffusion and its role in active deformable media

TL;DR

This paper introduces a new model coupling diffusion tensors to mechanical stress in active deformable media, demonstrating how stress influences anisotropic diffusion patterns and wave propagation, with potential applications in bio-materials like the heart.

Contribution

It presents a novel theoretical framework and experimental validation for stress-driven anisotropic diffusion, advancing understanding of mechano-electrical feedback in deformable biological tissues.

Findings

Stress couples to diffusion tensors, creating anisotropy.

Numerical simulations show effects on wave conduction and drift.

Experimental results support the model's predictions.

Abstract

We propose a new model to describe diffusion processes within active deformable media. Our general theoretical framework is based on physical and mathematical considerations, and it suggests to use diffusion tensors directly coupled to mechanical stress. A proof-of-concept experiment and the proposed generalised reaction-diffusion-mechanics model reveal that initially isotropic and homogeneous diffusion tensors turn into inhomogeneous and anisotropic quantities due to the intrinsic structure of the nonlinear coupling. We study the physical properties leading to these effects, and investigate mathematical conditions for its occurrence. Together, the experiment, the model, and the numerical results obtained using a mixed-primal finite element method, clearly support relevant consequences of stress-assisted diffusion into anisotropy patterns, drifting, and conduction velocity of the resulting excitation waves. Our findings also indicate the applicability of this novel approach in the description of mechano-electrical feedback in actively deforming bio-materials such as the heart.

Figures

Tables

Table 1. Table 1: Model parameters specifying the phenomenological RD cardiac dynamics and the active-stress mechanics ( 2.1 )-( 2.5 ).

$\Omega ={(0,250)}^2$, $a=0.1$, $b=0.15$, $k=8$, $ϵ=0.001$, ${ϵ}_0=0.1$,

${\mu }_1=0.12$, ${\mu }_2=0.3$, $k_{Ta}=9.58$, $c_1=6$, $c_2=2$, $D_0=1$,

$V_{\text{init}}=0$, $r_{\text{init}}=0.4$, ${\left(T_a\right)}_{\text{init}}=0.2$, $D_1,D_2$ from Figure 2

Equations

$$\frac{\partial V}{\partial t}$$

$$\frac{\partial T_a}{\partial t}$$

$$I_{\mathrm{ion}}$$

$$f(V,r)$$

$${\sigma }_{ij}$$

$$d_{ij}({\sigma }_{ij})$$

$$[{\sigma }_{ij}]=\left[\begin{array}{cc}{\sigma }_1& 0\\ 0& {\sigma }_2\end{array}\right],[{\sigma }_{ik}{\sigma }_{kj}]=\left[\begin{array}{cc}{\sigma }_1^2& 0\\ 0& {\sigma }_2^2\end{array}\right],$$

$$[d_{ij}]=\left[\begin{array}{cc}1+D_1{\sigma }_1+D_2{\sigma }_1^2& 0\\ 0& 1+D_1{\sigma }_2+D_2{\sigma }_2^2\end{array}\right],$$

$$\frac{-D_1\pm \sqrt{D_1^2-4D_2}}{2D_2}\ge 0,\text{and}{D}_1\ge \pm 2\sqrt{|D_2|}.$$

$$\mathbf{u}={\left[50{\sin }^2\left(\frac{\pi }{50}[t-t^{\ast }]\right),0\right]}^T,$$

Full text

Key words: Active deformable media, Stress-assisted diffusion, Reaction-Diffusion, Electro-Mechanics, Finite elasticity, Cardiac dynamics.

PACS numbers: 82.40.Ck, 77.65.Ly, 87.19.Hh, 87.10.Pq.

1 Introduction

Excitable media, whether of biological type or not, represent complex nonlinear systems. They are often of electrochemical nature, and can typically be coupled to several multi-physical factors as heat transfer or solid and/or fluid mechanics. A remarkable example is the heart, where nonlinear bioelectrical waves propagate on a complex anatomical background undergoing large mechanical deformations and facing strong interactions with biological fluids. More precisely, cardiac contraction results from the combination of a complex emerging behaviour where subcellular ion dynamics induce the overlapping of protein filaments, rapidly scaling up to both the cellular and tissue scales through a process known as excitation-contraction coupling and, as main topic of the present work, its reverse effect known as the mechano-electric feedback (MEF) [20, 22]. Studying the spatiotemporal dynamics of excitation waves in the heart is of paramount importance in the understanding of a large class of processes including depolarisation, repolarisation and period doubling bifurcations occurring in the transition towards chaotic regimes (arrhythmias) [50, 11, 6].

Still in the context of cardiac dynamics, a large number of experimental data is available to describe ionic and electrophysiological processes at many spatio-temporal scales. However, due to diverse technical reasons, a common practice is to biochemically suppress any mechanical feedback to record these data, implying that any back-reacting effect intrinsically due to electromechanical interactions is systematically neglected. Nevertheless, the importance of understanding the interplay between the reaction-diffusion (RD) dynamics with mechanical deformation is quite clear, and a recent growing interest in refining a companion mathematical model for the dynamics of higher complexity models has been observed [37, 38, 41, 42]. Even though several subcellular contributors to cardiac MEF have been extensively studied (as for instance, stretch-activated ion channels), their proper and consistent integration into tissue-level models has remained a challenging task; further having a very limited clinically-translatable application and validation.

In the last fifteen years, the relation between cell or tissue-based electrophysiology models with mechanical deformation of soft tissues has been formulated in terms of active stress [30], active strain [9, 32, 15], with the addition of stretch-activated currents [35, 36, 39, 51], or recently combined with membrane capacitance changes [12]. However, in all of these works the key physical ingredient ruling the spatio-temporal dynamics of the membrane potential, i.e. the diffusion tensor (in this context, the conductivity), has been systematically considered independent of mechanical deformation, even if general constitutive prescriptions suggesting a possible interaction were advanced [30]. Processes related to MEF have a fundamental role in a wide variety of passive physical systems. Notable examples comprise corrosion [56], rock anisotropy [19], glass transition [10], dissolution phenomena [28], electromigration [52], hydrogen trapping [7], as well as swelling effects [16]. Clear evidence for the existence of such a coupling in biological systems has also been recently observed in strain-dependent oxygen diffusivity in cartilage [18, 55], and in transcription factors within the cell nucleus [31]. Regarding the specific context of active biological media, connections forming gap junctions in cardiomyocytes (and associated to intercellular communication and mesoscale diffusion) have been recently discussed in terms of their mechano-sensitive properties [46]. Furthermore, a quantitative analysis on the specific effects of stretch into connexins in terms of hemichannels has been experimentally verified in a number of different cellular preparations [3, 8].

In the perspective of the present work, an important number of experimental studies have demonstrated key MEF effects in ventricular myocardium and atrial tissue (see e.g. [42, 41] for an extended review). Specific applications include atria arrhythmias, where basically three main results are available. First, the spatiotemporal distribution of atrial excitation depends strongly on the anatomical substrate [43]. Secondly, upon mechanical loading (stretching), the conduction velocity of the excitation wave decreases, a beat-to-beat interval variability appears, early afterdepolarisations, ectopic excitations and a higher vulnerability to atrial fibrillation are present [25]. Third, atrial tissue undergoes multiple high frequency and unstable rotors (spiral waves) when subject to constant and variable stretch states [54]. Experimental evidence of MEF has also been studied during ventricular loading, indicating a strong relationship between variations in the conduction velocity and strain anisotropy [23, 14, 29, 12].

The fact that electrical properties of solids undergo intrinsic modifications due to (even infinitesimal) deformations has been a subject of study since a few decades (see e.g. the classical volume by Landau et al. on the electrodynamics of continuous media [24, pp. 69]). These changes suggest the representation of the corresponding models using strain-dependent (or also stress-dependent) dielectric tensors. In particular, this dependency in turn affects the electromagnetic dynamics by enforcing inhomogeneous and anisotropic patterns in structures that were not necessarily so. Experimental evidence supports the main present assumption that diffusion depends on deformation. Therefore, and thanks to first principles, one can readily Taylor-expand a given diffusion tensor in terms of the deformation quantities.

In this note we present a novel formulation for the description of soft active deformable media within the context of coupled reaction-diffusion-mechanics systems, and employ nonlinear cardiac dynamics as a main motivating example. At this point we stress that the concept of stress-assisted diffusion has been originally formulated for generalised composite media (see [1, 53, 21, 27] and the references therein), but many resemblances exist with respect to the active deformation of soft tissues. Most notably, here we have found that an anisotropic and inhomogeneous diffusivity is naturally induced by mechanical deformations, thus affecting the nonlinear dynamics of the spatiotemporal excitation wave. This important fact implies that the present formulation can recover and generalise a large number of electromechanical models based on basic FitzHugh-Nagumo-type descriptions [35, 2]. The most relevant additional parameters are here the weights accompanying the stress when incorporated into the diffusion tensors, and therefore we study the plausibility of specific choices in the model parameter space. Our assessment is conducted for stretched tissues, focusing on appropriate physical indicators as conduction velocity and propagation patterns, and also carefully identifying conditions leading to the stability of the coupled system.

Before stating the proposed set of governing equations describing the stress-assisted diffusion model, we provide a proof-of-concept, elementary experiment. It consists of a dye diffusing on a sponge, where stress-assisted diffusion is readily encountered by simply applying stretch to the body. Figure 1 presents two mechanical states of the system, where the effects of the underlying structural features of the medium (with respect to the expected propagation patterns without stretch) are illustrated. We do not claim that such a process faithfully represents the physical mechanisms ruling biological media, and especially not active materials. Nevertheless, some interesting features are indeed shared with the general framework we intend to analyse.

2 A stress-assisted electromechanical model

We centre our investigation on an active stress RD model describing non-oscillatory cardiac tissue that supports stable propagation of excitation waves [30, 35]. We frame our modelling into finite elasticity, where one identifies the relationship between material (reference) and spatial (deformed) coordinates, indicated by $X_{I}$ and $x_{j}$, respectively, via the smooth map $x_{i}(X_{I})$ that determines then the deformed position of a point $x_{i}$ originally located at $X_{I}$. We indicate with $J$ the Jacobian of the map. In the deformed configuration the proposed equations read:

[TABLE]

with constitutive prescriptions for the RD system [35]

[TABLE]

and for incompressible isotropic materials $(J=1)$

[TABLE]

Equation (2.1) provides a non-dimensional, two-variables RD model where $V$ here represents the transmembrane potential and $r$ is the recovery variable, whose dynamics is prescribed by Eqs. (2.3). Equation (2.2) provides the contractile (active stress) dynamics of the tissue and the balance of linear momentum (law of motion). Here $\epsilon(V)$ is a Heaviside function introducing a discontinuous coupling (a switch) with the dynamics of the voltage [30]: $\epsilon(V)=\epsilon_{0}$ if $V<0.05$, and $\epsilon(V)=10\epsilon_{0}$ if $V\geq 0.05$. Equation (2.4) provides the total equilibrium stress in the current deformed configuration ($\sigma_{ij}$ is the Cauchy stress [48]). Such a formulation highlights the two main contributions due to the multiscale behaviour of the tissue, i.e. (i) the passive stress component originating from the classical strain energy function adopted for incompressible Mooney-Rivlin materials [33] and (ii) the active stress component following the formulation proposed in [30], contributing as an additional hydrostatic term to the total stress tensor only. The second order tensors $\delta_{ij}$ and $b_{ij}$ are the identity and left Cauchy-Green strain [48], respectively. Model parameters are provided in Tab. 1.

Equation (2.5) introduces the stress-assisted diffusion equation [1] generalising the effects of tissue deformation on action potential spread for isotropic and incompressible nonlinear elastic solids. On these assumptions, and according to the representation formula for isotropic second order tensors [48], the diffusive flux accounts for hydrostatic stress effects on the diffusing species and is characterised by three contributions: $D_{0}\delta_{ij}$ recovering the classical linear diffusion equation [30], while $D_{1}\sigma_{ij}$ and $D_{2}\sigma_{ik}\sigma_{kj}$ allow for stress-induced anisotropies into the diffusion tensor (this part is associated to the theory by Landau et al. [24] for infinitesimal deformations). We disregard here stretch-activated currents, as introduced later in [35], and follow instead the treatment in [30] to highlight the new main ingredients of the work.

3 Numerical results

Numerical method.

The discretisation of (2.1)-(2.5) follows the mixed-primal finite element formulation proposed in [45], where strains and stresses are computed directly without postprocessing them from low-order discrete displacements. This approach is preferable in the present context as we aim at analysing properties of the stress and diffusion tensors directly. Consequently, stresses are approximated with row-wise lowest order Raviart-Thomas elements [44], displacements with stabilised Brezzi-Douglas-Marini elements of degree one [4], whereas pressure (if actually needed) can be recovered from the discrete stress. The remaining fields (transmembrane potential, recovery variable, and active tension) are approximated with continuous, piecewise linear elements. We employ unstructured triangular meshes of different resolutions, and we set a fixed time step $\Delta t=10^{-3}$. The solver consists of an outer time advancing loop (of first order backward Euler type), an embedded fixed point iteration to decouple the electrophysiology RD system from the finite elasticity discrete equations, and inner Newton steps for the solution of the nonlinear mechanics. In turn, the AP model (2.1) is solved semi-implicitly (taking only the reaction terms explicitly), determining a formal CFL condition associated to the finest mesh and supporting the choice of the timestep mentioned above. All linear solves are carried out with a PETSc built-in LU solver using ILU factorisation as preconditioner.

Parameter space.

In order to determine admissibility regions for the modified RD system, we conduct an inspection of the ellipticity regimes depending on the model parameters (see e.g. [47]). Let us first consider a uniaxial tensile configuration of a two-dimensional domain where we induce a target pattern (see Figure 3, left) and identify horizontal ($X_{0},X_{1}$) and vertical ($Y_{0},Y_{1}$) positions. Fixing the stretching of the domain at the maximum physiological level for cardiac tissue, i.e. 20% of the resting length, one can obtain the maximum tensile components and a diagonal form of the Cauchy stress tensor can be readily derived in the central uniform tensile region [33], i.e.:

[TABLE]

where $\sigma_{1},\sigma_{2}$ correspond to the eigenvalues of the stress tensor that overlap with the two diagonal components of the stress in the Cartesian coordinate system, i.e., $\sigma_{11},\sigma_{22}$. Accordingly, the generalised diffusion tensor (2.5) becomes

[TABLE]

where a second order polynomial expression is obtained and parametrised over the two coefficients $D_{1},D_{2}$. We can then readily identify the regimes where the modified conductivity is a positive definite tensor. This translates in requiring that the graph of $y=1+D_{1}\sigma+D_{2}\sigma^{2}$ is always positive, i.e.,

[TABLE]

In particular, this condition ensures that both tensile and compressive states are physically allowed according to the material response without violating ellipticity.

We proceed to characterise the theoretical parameter space for an idealised tensile state $(\sigma_{1}=\sigma_{2}=1)$. Figure 2 displays the map $f(D_{1},D_{2})=1+D_{1}/D_{2}+1/D_{2}$, highlighting in blue colour the admissible regions where ellipticity holds. On the other hand, we perform numerical simulations (with stresses $\sigma_{1},\sigma_{2}$) and extract the computed values of the stresses at selected locations ($X_{0},Y_{0},X_{1},Y_{1}$) whose ellipticity conditions are shown in Figure 2. Finally, the physical parameter space is derived directly from numerical simulations for the same uniaxial condition in Figure 2. We compute the conduction velocities of the generated wave along the selected locations, $CV_{x},CV_{y}$, and identify regions with loss of ellipticity when the travelling waves present a propagation fault (non coloured region). In particular, the produced values of $CV_{x}/CV_{y}$ are for the most larger than 1 since the mechanical stretching is applied in the $x-$direction and the computations are non-diverging for most positive parameter values. Close to the region $D_{1}=D_{2}=0$, the CV ratio is $\sim 1$, thus recovering the conduction velocity expected for an isotropic medium (that is, the same independently of the direction). Particular combinations with $D_{1}$ or $D_{2}<0$ are characterised by $CV_{x}/CV_{y}<1$ thus implying a reduction of the conduction velocity in the direction of stretching (see regions in blue, in Figure 2). In contrast, for larger positive parameter values the wave moving in the $x-$direction propagates much faster than the one in the $y-$direction (up to a fivefold in our analyses). In all the selected scenarios, the analysis clearly indicates a marked anisotropic effect on the diffusion process, due to stress.

Stress-induced anisotropy.

Figure 3 shows three representative examples of different levels of the stress-assisted anisotropic coupling. One can observe how the propagating behaviour is modified by the level of anisotropy induced on the tissue, initially homogeneous and isotropic, thanks to the stress-assisted mechanism. Based on the selected parameters, a faster propagation is obtained in the direction of loading where a horizontal uniaxial stress state is induced. In the first case (corresponding to Figure 3), the choice $D_{1}=D_{2}=0$ does not induce any effect on the excitation wave and the only visual difference is due to the geometric mapping from the reference to the deformed configuration of the computational domain. Actually, this setting (which is the expected behaviour in the classical active stress approach for incompressible isotropic materials [30]) is the most commonly employed in models and simulations for the electromechanical coupling.

The second case (see Figure 3) is characterised by $D_{1}=-0.01,D_{2}=0.01$ and shows a mild, though significative, anisotropic effect due to the applied mechanical stretching in the horizontal direction. The obtained behaviour is the one expected based on the few cardiac optical mapping data describing such a phenomenon in a moving frame [28] and, out of these three cases, it is qualitatively the most appropriate, in terms of physiological behaviour of cardiac dynamics.

The last run (see Figure 3) uses the values $D_{1}=-0.25,D_{2}=0.25$, and our results produce a $CV_{x}$ about a threefold higher than $CV_{y}$. These tests point out a critical example of the effect of mechanical stretching on the diffusive properties of the species in terms of transmembrane potential $V$ (top panels) and recovery variable $r$ (lower panels). The observed level of distortion induced on the excitation wave may be expected in some pathological scenarios with drastically compromised excitation-contraction dynamics [34]. We remark that this analysis is of particular importance in the study of the effects of the electromechanical coupling into spiral wave drifting phenomena and atria arrhythmias.

Stress-induced arrhythmogenesis.

Additional effects on arrhythmogenesis are exemplified with a simulation of sustained spiral dynamics. Figure 4 compares the long run behaviour of an equally induced spiral within a uniaxially stretched tissue for three increasing levels of the stress-assisted diffusion parameters $D_{1}$ and $D_{2}$. For the chosen parameters’ set from the original model a stationary circular meander pattern is expected (Figure 4), although spiral drifting is observed as the stress-assisted diffusion is activated (Figure 4). In addition, totally distorted excitation waves are observed when a strong anisotropic coupling is enforced (Figure 4) via higher values of the parameters $D_{1}$ and $D_{2}$. The corresponding recovery variable is also provided and confirms the highlighted differences.

The analysis is further quantified by rendering the modified diffusion tensors in terms of their tensorial glyph representation (see e.g. [13]). These plots (see Figure 5) are generated by solving local eigenvalue problems using the tensor at hand and displaying an ellipsoid whose shape and size determine the magnitude of eigenvalues and direction of eigenvectors and indicating the local patterns of anisotropy and inhomogeneity. For reference, we also depict the $L^{2}-$norm of the resulting diffusion tensor in each case, that is $D_{0}C_{IJ}^{-1}$, $D_{IJ}$, and $S_{IJ}$, respectively. As expected, constant diffusion, Figure 5, and stress-assisted diffusion, Figure 5, lead to different representations of the eigenproblem. Even if we only account for the uniaxial case (represented in Figure 5), substantial differences are expected when multiaxial stress patterns are considered.

According to [35], where stretch-activated currents are present, we observe a small spiral drift velocity of about 1.4% variation with respect to pulse wave velocity. This means that several spiral rotations are required to reach the boundary of the domain such to detect the exit of the excitation wave from the domain (auto-defibrillation). However such small differences can lead to substantially different scenarios at long times, as observed in strongly nonlinear dynamical systems  [5, 17]. This behaviour is in line with the expected vulnerability to atrial flutter and atrial fibrillation during atrial dilatation [43, 25, 54]. At this stage we remark that our theoretical study does not imply a stretch-activated current contribution (still debated in the literature [41]) but only relies on the effect induced by the generalised diffusivity tensor.

Enhanced effects due to mechanical pacing.

Finally we carry out a set of simulations of the fully coupled stress-assisted diffusion model, concentrating in the case where a periodic and uniaxial mechanical loading is applied on the right end of the square slab $\Omega=(0,250)^{2}$. Again, the material and electrical parameters are taken from Table 1. The boundary conditions for the nonlinear elasticity problem correspond to zero normal stresses on the horizontal boundaries, a clamped left boundary, and on the right we prescribe the displacement

[TABLE]

where $t^{*}=100$ indicates the onset time for the mechanical loading. We use a S1-S2 protocol to initiate spiral wave dynamics, and focus our attention on the long term behaviour of the patterns produced according to three different levels of stress-assisted diffusion. In Figure 6 we portray the trajectory of the spiral tip (also showing a zoom on the relevant region) over the time interval $t\in[400,1000]$. The location of the spiral tip is recorded by computing a piecewise constant field defined as the module of the gradient of the recovery variable, and then obtaining its global argmax over the spatial domain. Our results indicate a marked progressive drifting of the spiral tip, which increases as the relevant parameters depart from the stable values.

4 Discussion

In this work we propose a generalised MEF formulation for active deformable media based on the concept of stress-assisted diffusion. The leading mechanism of the modification of wave propagation corresponds to the induction of anisotropy in an originally isotropic medium. In this study, we simulated a single spiral wave that drifts or degenerates according to the amount of stress-assisted coupling. The proposed physics, experimentally and theoretically based on physical grounds, enriches the phenomenological description of electro-mechanical RD systems in the context of cardiac dynamics and excitable deformable media in general. The numerical results confirm that the generalised formulation induces anisotropy in the wave propagation starting from an isotropic medium. Secondly, the RD dynamics are modified, as clearly evidenced by the change in spiral meandering (i.e. drift of the spiral tip [30, 35]). In addition, the generalised model leads to inhomogeneous and unpredicted excitation waves, especially pronounced for extreme values of the coupling coefficients as observed in long run tests and under mechanical pacing.

To close this section we collect a few possible extensions of this work. The connection between spatio-temporal variation of the tissue anisotropy and the identification of model parameters in general cases remains unclear, but it is by no means a trivial task. One would definitely require to analyse deformation patterns interacting with a number of intrinsic properties (not considered here) as domain curvature, a priori anisotropy, pre-stretch, structural heterogeneities, and many others.

Also, we would be interested in analysing important microscopic effects related to gap junction proteins density at cell-cell interfaces. These properties have an important effect on the macroscopic diffusivity and its relation with the conduction velocity [49], especially under pathological remodelling conditions. Other envisaged modifications include different mechanical constraints [26], and the incorporation of beat-to-beat adaptation (see the recent review [40]).

An important complement to our generalised MEF formulation would be the characterisation of admissibility ranges for an elementary mechanical setup, together with the experimental validation of the proposed model. The latter could be carried out by fine-tuning stress-assisted diffusion parameters via optical mapping of cell cultures, myocardial tissue slices, and whole heart geometries (with particular attention to atrial tissue).

Several other applications are foreseen, as stretch-activated currents and MEF in cell-cell coupling phenomena. However a deeper understanding of the interacting regimes (also including model parameters and biophysical descriptions) is still required before engaging in the study of these mechanisms using homogenisation theory.

Acknowledgements.

This work has been partially supported by the Italian National Group of Mathematical Physics GNFM-INdAM, by the International Center for Relativistic Astrophysics Network ICRANet, by the London Mathematical Society through their Grant Scheme 4, and by the EPSRC through the Research Grant EP/R00207X/1.