Finite element methods for large-strain poroelasticity/chemotaxis models simulating the formation of myocardial oedema

In this paper we propose a novel coupled poroelasticity-diffusion model for the formation of extracellular oedema and infectious myocarditis valid in large deformations, manifested as an interaction between interstitial flow and the immune-driven dynamics between leukocytes and pathogens. The governing partial differential equations are formulated in terms of skeleton displacement, fluid pressure, Lagrangian porosity, and the concentrations of pathogens and leukocytes. A five-field finite element scheme is proposed for the numerical approximation of the problem, and we provide the stability analysis for a simplified system emanating from linearisation. We also discuss the construction of an adequate, Schur complement based, nested preconditioner. The produced computational tests exemplify the properties of the new model and of the finite element schemes.

of inflammatory oedema in the context of blood-brain barrier failure [50], and immune systems for small intestine [75,82] as well as for myocarditis [39].
In this work we are concerned with the latter application involving the formation and evolution of oedema, a build up of excess of fluid content in the intercellular space, in myocardial tissue. Local infection of tissue by pathogens contribute to an inflammatory reaction driven by the immune system with the aim of protecting the compromised region against invading organisms. Driven by oncotic and hydrostatic pressure gradients, there is a transvascular water flux: capillaries transport fluid into and out from the interstitium (the intercellular region), and lymphatic nodes contribute to returning fluid out from that space. This work also considers that, before infection occurs, the body is under osmotic equilibrium since the capillary barrier is permeable to ions. In this case, ions have no significant effect on the interstitial flow under normal conditions, otherwise, a previous disease may exist. The presence of pathogens inside a cell as well as tissue damage caused by them may change this scenario, i.e., ion-concentration may change, impacting local pressure. However, in this work, we assume, for simplification purposes, that the main cause of oedema formation is the presence of the pathogen in the tissue and the immune response to it. In the process of inflammatory response, local vasodilation occurs, which not only allows leukocytes to leave the bloodstream and to access the infection site but also increases the accumulation of fluid in the intercellular space, see Figure 2.1. However such abnormal accumulation of interstitial fluid content may lead to local stress generation, tissue deformation and swelling, eventually producing a number of complications that can impair the normal function of the tissue. Oedema may be found in several types of tissue. Heart oedemas are commonly observed in myocarditis and they constitute one of the main criteria used in medical practice to its diagnosis [41,62]. Deriving refined continuum-based models for cardiac oedema formation valid in the large strain regime and that include 3D personalised geometries have the potential to describe mechanistic processes that are observable by clinical evidence and to provide deeper insight into the complex multiscale mechanisms that are inherent to the impairing of the healthy cardiovascular function.
We extend the preliminary results from [39,40] and develop a simplified phenomenological model for the dynamic interaction between poroelastic finite strain deformations and the chemotaxis of leukocytes towards pathogens. The present framework assumes that the porous skeleton admits heterogeneities at the macro scale, and that the pores are all interconnected. We also consider that inertial forces are negligible. The solid-fluid interaction describing the distribution of flow within the deformable skeleton is then cast in Lagrangian coordinates using the solid displacement, the fluid pressure, and the nominal porosity (ratio between the pore volume and the total volume), whereas the immune system dynamics are represented by the concentrations of leukocytes and a pathogen. The model for the immune system has also been modified to better represent the dynamics that occur after the pathogen has been eliminated. Furthermore, the mechanical feedback into pathogens and leukocyte populations are now taken into account, and we have explored these aspects further in our recent work [52], which uses the present formulation for the coupled poromechanical/chemotaxis equations.
Even if many numerical methods for the approximation of linear poroelasticity and their convergence analysis can be found in the recent literature (see [12,[17][18][19]55] and the references therein), much less attention has been given to formulations addressing large-strain poromechanics. In this regard, recent works include an enriched Galerkin framework [27], stabilised finite elements [15,84] and hybrid finite elements as well [83]. Although most contributions address the problem in a Lagrangian frame of reference, some alternative formulations include descriptions in Eulerian [65] and ALE [22] coordinates. Besides the stress response, nonlinear models differ from the linear ones in that pressure is not a primary variable, but instead it is given by constitutive modelling [25,31]. Instead, we follow the framework of [30,53] (and other similar models) where the issue is circumvented through solid phase incompressibility, which yields a Lagrange multiplier acting as hydrostatic pressure.
The extension of such solution techniques for large scale simulations is highly non-trivial due to the requirement of efficient preconditioned iterative methods for the poroelasticity problem. Preconditioning strategies are mainly given by operator preconditioning [9,46] and block factorisations [42,79]. In this regard, block preconditioners are more easily generalisable, but an efficient implementation requires an adequate approximation of the Schur complement, usually approximated through a fixed-stress formulation in poroelasticity [38,42,79]. Despite all of this, an efficient numerical solver for a linearisation of large-strain poromechanics with nominal porosity as a primary variable remains largely unaddressed.
In the present case we formulate the set of nonlinear equations using as primary variables the solid displacement, the fluid pressure, the Lagrangian porosity, and the concentrations of leukocytes and pathogens. A finite element method is employed for the space discretisation [67], combined with a first-order backward Euler time advancing scheme. The method uses a MINI element for the displacement-porosity pair [6] (see also [8,23] for the case of hyperelasticity), and linear Lagrangian elements for the remaining unknowns. As the stability of the fully nonlinear system is still a paramount task (and an open problem even for the segregated nonlinear poroelasticity without the reaction-diffusion coupling, see, e.g., [12]), we address the solvability of the linearised set of equations, focusing on the semidiscrete in-time formulation. For the corresponding analysis, we propose a fixed-point strategy, and then, Schauder fixed point theorem, together with the Fredholm alternative, the Babuška-Brezzi theory, the Lax-Milgram lemma, and classical theory of quasi-linear equations, allow us to assert unique solvability robustly with respect to the material parameters of the linearised poroelastic solid. Then, we address the numerical approximation of the linearised problem through a Krylov subspace method with a Schur complement preconditioner.
The main advantages of the proposed mathematical model and the associated computational methods are i) a consistent theoretical framework valid for finite strains that stems as the natural generalisation of three-field linear poroelasticity from [60] (which is written in terms of displacement, fluid pressure, and total pressure), ii) the versatility of the formulation to accommodate 2D or 3D geometries, iii) the accuracy of the numerical scheme and efficiency of the Krylov subspace iterative methods under consideration, and iv) the potential of the methodology in replacing current invasive methods for the detection of interstitial fibrosis and myocarditis (such as endo-myocardial biopsy) by techniques hinging only on MRI data.
We have organised the contents of this paper in the following manner. Section 2 outlines the main details of the model problem, motivating each component in the balance equations and stating main assumptions together with typical boundary conditions. In Section 3 we present the weak formulation of the nonlinear coupled problem, as well as the derivation of the linearised form of the weak formulation. Then, in Section 4, we address the solvability of the linearised set of equations. Section 5 contains the derivation of a finite element scheme, including the fully discrete problem in matrix form, and we address the construction of efficient solvers tailored for the multiphysics coupled system. In Section 6 we present a number of computational results, consisting in a simple numerical study concerning the sensitivity analysis of relevant model parameters, the verification of spatio-temporal convergence, and an investigation of different cases on simplified and physiologically accurate geometries. We close with some remarks and a discussion on model extensions collected in Section 7.

Continuum model and proposed set of field equations
The different scales and structural components in the tissue suggest to employ a continuum model. We assume that the tissue is a poroelastic medium with interstitial fluid completely filling the void spaces. Adopting usual kinematic conventions and notation, let us consider a domain Ω ⊂ R d , d = 2, 3 representing the volume occupied by a deformable porous structure in its reference configuration. We will denote by n the outward unit normal vector on the boundary ∂Ω. We also disjointly partition the boundary into the subboundaries Γ and Σ (also in the reference undeformed state) where different types of boundary conditions, associated with the equations of motion, will be applied. A material particle in the undeformed domain Ω is identified by its position x (and, unless specified otherwise, all differential operators such as gradients and divergences will be evaluated with respect to these undeformed coordinates), whereas u : Ω → R d will denote the displacement field defining its new position x + u in the deformed configuration. The tensor F := I + ∇u is the gradient (applied with respect to the fixed material coordinates) of the deformation map and I denotes the identity tensor; its Jacobian determinant, denoted by J = det F = det(I + ∇u), measures the dilation (solid volume change during the deformation); and C = F t F is the right Cauchy-Green deformation tensor, where the superscript () t denotes the transpose operator.
The non-viscous filtration flow through the non-deformed porous skeleton can be described by Darcy's law [2,3]. Adopting its parabolic formulation solely in terms of pressure head (which condenses momentum and mass conservation equations for the fluid) and combining it with the conservation of linear momentum in the solid, we obtain the extension of Biot consolidation equations to the regime of large strains. This can be carried out using the descriptions already available in, e.g., [15,53,84]. In this paper we choose to write them using the first Piola-Kirchhoff stress tensor P = P eff − αpJF −t and the effective (hyperelastic) stress P eff = ∂Ψs ∂F . For incompressible neo-Hookean and Holzapfel-Ogden solids one has the following forms (see, e.g., [45,63]) , s, f s} are material parameters and we have used the notation (ζ) + := ζ if ζ > 0 or zero otherwise, and the invariants and pseudo-invariants of C measuring isotropic and direction-specific stretches I 1 (C) = tr C, Here f 0 and s 0 denote local alignment of fibres and sheetlets which is inherent to the local microstructure of the myocardial tissue [58]. Histological studies based on diffusion tensor imaging indicate a fibre architecture that rotates 120 • from the epicardial to the endocardial surface. A volumetric part of the energy is added for the case of nearly incompressible materials, for example, the term U (J) = 1 2 λ s (lnJ) 2 (with λ s another material parameter) which is designed to enforce convexity or polyconvexity of the strain density [85].
Provided that the volumetric stiffness of the solid is substantially smaller than that of the constituent (which is reasonable for soft living tissues, see, e.g., [50,54]), volume changes of solid constituent are negligible, implying that where φ f − φ 0 is the nominal porosity change in the reference configuration [53] (here we are assuming that the mixture is saturated, that is, φ f + φ s = 1, where φ s is the volume fraction of the solid). Recall that the Lagrangian porosity is the natural state variable being work-conjugate to the fluid pressure p (see, e.g., [59]). Should pre-stress be considered, we need to take instead J = J 0 + φ f − φ 0 , where J 0 is the volume associated with the pre-stress. Note that the specific form of the strain energy density Ψ s is considered irrespective of the fluid saturating the solid. It can also be derived from subtracting the volumetric free energy of the fluid from the total Helmholtz free energy of the mixture Ψ s = Ψ −φ f Ψ f , where φ f is the nominal (Lagrangian) porosity measuring the current fluid volume per unit reference total volume [59,84]. Note also that rearrangement of tissue components as a result of change in fluid content will imply an additional solid deformation as well as a stress modification [27,50,61].
The problem is then stated in a Lagrangian framework and in terms of the solid displacement u, the fluid pressure p, and the nominal porosity φ f , as where D t stands for the material time derivative, and Div and Div are the divergence operators (with respect to the coordinates in the reference configuration) for vector and tensor fields, respectively. Equation (2.2a) is the balance of linear momentum for the multiphase system expressed in terms of the true total stress; (2.2b) states the mass conservation (or continuity) of fluid content in a material framework, written in terms of the nominal flux that includes the movement of the fluid phase due to pressure gradient under Darcy's law (which in particular discards shear flow effects and is valid for small Reynolds numbers); and (2.2c) represents the aforementioned kinematic constraint relating volume and porosity. It states the material incompressibility of the constituents, but it still allows for compression of the medium since the pore volume can change locally in connection with modifications in the macroscopic mass density of the poroelastic material. A more complete discussion can be found in [53]. This description is independent of the material constitutive model for the solid matrix, and the balance of angular momentum simply furnishes the condition FP t = PF t . Here D t indicates the material time derivative (with respect to skeleton particles), b is a vector of external body loads, is a distributed source/sink strength term (local source volume of fluid added per volume of tissue per unit time) that accounts for fluid being added or removed from the interstitium by the blood and lymph capillaries and its specific form will be made precise later on, ρ s is the reference density of the porous matrix, ρ f is the reference density of the interstitial fluid, α is the Biot-Willis modulus assuming values between 0 and 1, and µ f is the dynamic viscosity of the interstitial liquid. In the inflammatory response process, local increment vascular permeability occurs, which allows leukocytes to leave the bloodstream to access the infection site which also leads to fluid accumulation in the interstitial space. Driven by oncotic and hydrostatic pressure gradients, blood capillaries transport fluid into and out of the interstitium, and lymphatic capillaries contribute to the resolution of the oedema.
The porosity-dependent heterogeneous tensor of intrinsic permeability (or hydraulic conductivity) κ can be defined by using normalised Kozeny-Carman (KC), exponential (E), or by power-law (PL) relations defined, respectively, as is the typical diameter of pores or of solid grains, and M 0 is the slope of the normalised permeability curve [7,14,84]. As porosity is a volumetric quantity, these forms of permeability are isotropic and therefore indifferent to rotation and shear effects. A relatively simple phenomenological extension that accommodates anisotropy due to tissue microstructure is given as follows (using as example the Kozeny-Carman law) where κ f 0 , κ s0 are weights for the principal hydraulic conductivities such that the ratio κ f 0 /κ s0 ≈ 2.5 coincides with the anisotropy ratio typically used for the conductivity tensor in the extracellular domain of bidomain models [29,74] (and much milder than the ratio observed in the intracellular domain, or in other porous skeletons such as sand or cartilage [21,37]). This allows to reproduce the preferential perfusion direction of the micro-vascular system following the orientation of cardiac fibres f 0 , s 0 . Other transformations that go beyond the kinematic nonlinearity induced by change of frame can be found in models exploiting strain-dependence [7] or stress-assistance mechanisms [26]. Unless otherwise specified, we will use (2.3).
The infection process starts when an invader (pathogen) enters the poroelastic tissue. Once there, it starts to propagate. The exact mechanism adopted for increasing the pathogen population (replication, in the case of viruses, or reproduction, otherwise) is not essential to our model. In our simplification, γ p represents the rate at which the pathogen propagates. The innate immune cells (leukocytes) use their receptors to sense the presence of non-self molecules characteristic of pathogens. After the presence of pathogens is identified, leukocytes start to produce pro-inflammatory cytokines, which increase vascular permeability and recruit more leukocytes to the site of infection. The model does not include cytokines because the concentration of pathogens and leukocytes can indirectly represent their effect on vascular permeability and recruitment of leukocytes. Additionally, we assumed a simplified reaction term for the leukocytes since among these cells some have a very long life span when compared to the timescale of the simulations considered in this work. For example, macrophages can live from months to years [76]. Evidently, leukocyte death from apoptosis and pathogen phagocytosis create a sink that was not considered in the current work. Leukocytes adhere to the endothelial cells that line the vessels, a phenomenon called margination (Figure 2.1). The increase in vascular permeability facilitates leukocyte extravasation to the tissue. The model considers that leukocytes leave the bloodstream to the tissue with an extravasation rate equals to λ pl . The leukocytes move towards the site of infection due to the presence of chemoattractants. The chemotactic rate is given by χ cl . Finally, leukocytes can start the phagocytosis of pathogens, eliminating them at a rate equal to λ lp (Figure 2.1). We denote by c p and c l the concentrations of pathogens and leukocytes, respectively. Both populations can also diffuse in the interstitial fluid. Their dynamics, in the current configuration and with derivatives in terms of the deformed coordinates, are governed by where D l , D p are positive-definite diffusion tensors for the species within the interstitial fluid, and r l , r p are reaction terms adopting the following specific forms (see, e.g., [39]) withγ p ,λ lp , andλ pl being constant coefficients (pathogen growth rate, leukocyte phagocytosis rate, and leukocyte migration rate, respectively) [39]. Note that the advecting velocity corresponds to the filtration velocity (that between the solid and fluid velocities, as the species are supposed to be transported in the interstitial fluid). Pulling back these advection-reaction-diffusion equations to the reference frame we obtain a modified system with kinematic nonlinearity Note that, in order to make the advection completely absorbed by the total derivative we have considered the total fluid-solid velocity. In turn, this implies that we are assuming that pathogens and leukocytes are able to move in the solid also.
Next we return to the continuity equation (2.2b) to specify the source (rate of fluid movement) . This term will encode the immune response feedback into the hydromechanics by setting a Starling-Hill function of the type The remaining coefficients 0 , p c , π c , π i , v max , k m , n, p 0 , S/V, L p0 , L bp , L br (normal lymph flow, capillary basal pressure, capillary oncotic pressure, interstitial oncotic pressure due to plasma protein, maximal lymph flow, half-live of pressure to lymph flow, Hill coefficient, normal interstitial fluid pressure, vessel area per volume unit, healthy tissue hydraulic permeability of the capillary wall, intensity of the pathogen infection into permeability, and influence of pathogens in the reflection coefficient; respectively) are positive model constants [39].
To close the system composed by (2.2a)-(2.2c) and(2.4a)-(2.4b), we need to provide suitable initial data and boundary conditions. We suppose that the system is initially at rest and with an homogeneous distribution of the leukocytes and pathogens in the domain In addition, on the boundaries Γ and Σ we prescribe traction and displacement, respectively; and zero flux everywhere on the boundary for the species interacting in the immune system. That is, we consider the following set of boundary conditions where t Σ is a given boundary traction applied in the current configuration.

Weak formulation and consistent linearisation
We derive a weak form for the coupled system, carry out a Newton-Raphson linearisation, and then discuss the solvability of a simplified system arising from the linearisation with respect to conditions at rest. That resulting problem coincides in structure with the coupled Biot/reaction-diffusion equations.

General weak form and in-time discretisation
In view of the essential boundary conditions for displacement and fluid pressure (2.6a),(2.6b), we define the Hilbert spaces associated with the classical norms in H 1 (Ω), and H 1 (Ω), respectively. The nonlinear weak form of the equations of motion and of the transport for the immune system can be established following a standard approach, that is, multiplying each field equation by a suitable test function, integrating over Ω, and invoking the divergence theorem whenever appropriate. To lighten the presentation, we drop the measures dx and ds from the integrals, which are to be understood throughout the manuscript as the Lebesgue measure. This leads to the following continuous, five-field weak formulation: For almost all t > 0, find (c p , We do not address time regularity in this work, and simply assume that the solution is weakly differentiable and continuous in time, the latter in order to obtain well-defined initial conditions. Next we partition the interval [0, t final ] into N evenly spaced non-overlapping sub-intervals of fixed length ∆t and apply a time semi-discretisation of (3.1) using the unconditionally stable, backward Euler's method. Starting from initial data, at each time iteration n = 0, . . . , N − 1 we have a nonlinear problem with solution (c n+1 , and φ n+1 f ∈ L 2 (Ω). Then we perform a Newton-Raphson method linearisation. This means that at each time step n + 1 we start the Newton iterates using as initial guess the solution at the previous time step c k=0 p = c n p , c k=0 l = c n l , u k=0 = u n , p k=0 = p n , φ k=0 f = φ n f , and for k = 0, 1, . . . we look for the solution increments (δc k+1 ,Ω , R k 3,Σ , R k 5 are vector and scalar residuals coming from the linearisation of (2.2a) and (2.2c) and containing terms associated to the previous Newton step k as well as body and traction loads. The terms R k,n 1 , R k,n 2 , R k,4 4 are scalar residuals from the linearisation of (2.4a), (2.4b), and (2.2b), respectively; which in addition contain contributions from the previous time step n.

Definition and preliminaries
For this section we draw inspiration from the stability analysis of a system similar to (2.2), recently proved in [12]; and from the study of the linearised hyperelasticity equations, recently studied in [36]. We restrict the description to Neo-Hookean poroelastic solids with material parameters µ s , λ s , for which the effective first Piola-Kirchhoff stress tensor is Note that a simplified version of the time-discrete tangent system (3.2) is readily obtained by focusing on the first Newton-Raphson iterate when starting from the following initial guess for displacement, fluid pressure, nominal porosity and species concentrations, which gives D t φ k f = 0, F k = I, J k = 1, ∇p k = 0. In addition, the residuals can be conveniently rearranged and we rescale the system in an appropriate manner using the parameters ρ f , ∆t, α, µ s , φ 0 , and λ s . Next we introduce the total pressure and the system can be recast in terms of u, ζ, p, c p , c l (that is, using ζ instead of φ f ), as follows where for a function scalar X n+1 , the expression δ t X n+1 stands for X n+1 −X n ∆t , and where we have renamed the bilinear forms, which are now defined as Here we have also set t Σ = 0. Note that the linear variational problem (4.2) corresponds to the linear and coupled reaction-diffusion / three-field Biot equations where the R i 's andR i 's are scalar and vector residuals, and the β i > 0 are additional constants arising from the linearisation of the reaction terms and from the change of variables (4.1). A variant of system (4.4) including advection and nonlinear reaction, has been studied in [77]. We follow that reference and in the remainder of this section we present the solvability, stability, and convergence analysis for the in-time problem associated with (4.2).
We will assume that the anisotropic permeability κ and the diffusion matrices D p , D l are uniformly bounded and positive definite in Ω. The latter means that, there exist positive constants κ 1 , κ 2 , and D min We also recall the following version of Korn's inequality [20]). This bound, together with the assumptions above on the permeability and diffusivity, imply the following coercivity and positivity properties . Moreover, thanks to the structure of the formulation, the bilinear form b 1 (which coincides with the usual non-diagonal bilinear form in the Stokes equations) satisfies the following inf-sup condition (see, e.g., [44]): For every ζ ∈ L 2 (Ω), there exists β > 0 such that Finally, we recall an important discrete-in-time identity and introduce the discrete-in-time norm respectively, which will be useful for the subsequent analysis.

Fixed-point scheme and unique solvability of the decoupled problems
With the aim to facilitate the comprehension and clarity of the forthcoming analysis, through this section and Sections 4.3 -4.5, we take up again the notation (·) n+1 given in Section 3.1 to denote the current time step of the involved variables. Moreover, in order to apply a fixed-point scheme, we modify the variational formulation (4.2) in an equivalent one: find (c n+1 respectively.
We then start with the fixed-point strategy. Let us define the operator S : where p n+1 ∈ H 1 Σ (Ω) and ζ n+1 ∈ L 2 (Ω) are the second and third components of the solution of problem: for a given c n+1 p ∈ H 1 (Ω). In turn, letS : for a given pair (p n+1 , ζ n+1 ) ∈ H 1 Σ (Ω) × L 2 (Ω). Therefore, we define the map T : and one readily realises that solving (4.7a) -(4.7e) is equivalent to seeking a fixed point of the solution operator T.
Now, in order to prove that the operator T is well-defined, we first need to investigate that the uncoupled problems defined by the operators S andS are well-posed. The fact that the poroelastic problem is well-posed for a given c n+1 p is given next. Its proof can be deduced by employing the Fredholm alternative approach and classical tools frequently used for showing the well-posedness of elliptic/parabolic equations. We refer to [  (4.14) In turn, the well-posedness of the uncoupled diffusion-reaction problem is given next. Proof. For each n, the unique solvability of (4.12a) -(4.12b) can be obtained by using the coercivity of the diffusivities D i i ∈ {p, l}, and classical tools for time-discrete approximation schemes of parabolic equations (see for instance [64]). Now, for the stability estimation, we begin by obtaining the corresponding result for  The same result can be derived for c n l . Finally, the estimate (4.15) follows by adding the aforementioned result to (4.16) and then, applying Gronwall's inequality, and using the resulting approach for n + 1.

Existence of a weak solution to the coupled linearised system
Let us define a closed subset of the Banach space L 2 (Ω) as for a given r > 0. We then restrict the space in (4.13) to the ball (4.17), and therefore, we look for a c n+1 In what follows, we prove the existence of that c n+1 p by means of the Schauder fixed-point approach. We start with a previous result. where C andC are positive constants depending on ∆t, µ f , µ s , κ 1 , β, α, with Proof. The proof starts by multiplying the equations (4.10a), (4.10b) and (4.10c) by 1 ∆t , α, and − 1 ∆t , respectively. Then, the estimate (4.18) follows after adding these resulting equations with v := u n+1 , q := p n+1 and ψ := ζ n+1 , respectively, using Cauchy-Schwarz and Young inequalities. We omit further details. Now, we are in a position to establish the existence of a fixed point of the operator T. This result is abridged in the following two lemmas.
for a given pair (p n+1 , ζ n+1 ) ∈ H 1 Σ (Ω) × L 2 (Ω), from which, we straightforwardly obtain  Finally, the desired result follows after substituting the estimate (4.18) into (4.21), and then, applying the assumption on data given at the statement of the present lemma.
We remark that there is no inconvenience with the second assumption on data given in Lemma 4.4 since it depends on a constant r that can be appropriately chosen (cf. (4.17)). Moreover, the first smallness assumption given in Lemma 4.4 is merely a theoretical condition to guarantee solvability of the continuous problem. However, even though for the constantsC andĈ defined by (4.19) and (4.22), respectively, the condition might not always be satisfied (especially for extreme values of the model constants), the numerical experiments obtained in Section 6 show satisfactory results even without the theoretical constraints. The present analysis could be improved by using a different approach that circumvents the current restrictions on data, but we are not aware of such a strategy being applicable directly to this context. Proof. We notice from the previous lemma, that T(K) is bounded in H 1 (Ω). In this way, the compact embedding of H 1 (Ω) into L 2 (Ω) together with boundedness of T(K) conclude that T(K) is relatively compact in L 2 (Ω). On the other hand, for the continuity property, we let {c n+1 p,k } k ∈ K be a sequence such that c n+1 p,k → c n+1 p in L 2 (Ω) as k → ∞, and as is known from (4.13), c n+1 p,k = T(c n+1 p,k ). In this way, proceeding as in [77, Lemma 2.8], we obtain that c n+1 p,k → T(c n+1 p ) in L 2 (Ω) as k → ∞, concluding the proof.
Finally, the following lemma concerning the existence of solution for problem (4.7a)-(4.7e), is merely an application of Lemmas 4.4 and 4.5, and the Schauder fixed-point theorem.

Stability of the linearised coupled problem
The following lemma establishes the continuous dependence on data for problem (4.7a) -(4.7e). Its proof follows similar arguments to those used in Lemmas 4.1, 4.2 and 4.7, and therefore we omit it here. -Although showing the solvability analysis of the linearised fully-continuous problem goes beyond the scope of this work, it is possible to extend the ideas presented here to establish that well-posedness using appropriate choices of Sobolev spaces, and performing a passage to the limit adapting to our problem the results from, e.g., [5] (which focus on reaction-diffusion-Brinkman problems). -The solvability analysis of the fully-discrete problem associated with (4.7a)-(4.7e) can be established similarly to the semidiscrete-in-time case. More precisely, with the finite elements spaces specified in (5.1), we can define an appropriate fixed-point operator and prove that it is well-defined thanks to the well-posedness of each uncoupled problem (see a general form in, e.g., [64]). Then, by employing the continuous dependence in the fully-discrete case (which follows exactly as in Lemma 4.8) in combination with the analysis from [5, Section 5.3], we can prove the continuity of the fixed-point operator described above. Finally, the desired result follows from an application of Brouwer's fixed-point theorem. -Given the approximation properties of the finite element spaces specified in (5.1) (and recalled in, e.g., [60]), and following the steps given by [77,Section 4], it is possible to derive an asymptotic O(h) convergence for the proposed method. This is also verified for the nonlinear case, as shown numerically in Section 6. -The convergence properties are also robust with respect to incompressibility of the solid phase.
-We stress that the arguments presented in this section are unfortunately not readily generalised to arbitrary hyperelastic materials, nor a generic initial guess for the Newton-Raphson loop. In all of the numerical tests performed below, this has not shown to be a problem in practice. Still, the lack of a global existence result makes it is possible for the model to be increasingly difficult in the simulation of other scenarios of interest, such as pathological ones.

Discretisation of the nonlinear problem
In order to define a Galerkin finite element method we denote by {T h } h>0 a shape-regular family of partitions ofΩ, conformed by tetrahedra (or triangles in 2D) K of diameter h K , with mesh size h := max{h K : K ∈ T h }. Given an integer k ≥ 1 and a subset S of R d , d = 2, 3, by P k (S) we will denote the space of polynomial functions defined locally in S and being of total degree up to k. Let us also denote by b K := ϕ 1 ϕ 2 ϕ 3 a P 3 bubble function in K, where ϕ 1 , ϕ 2 , ϕ 3 are the barycentric coordinates of the triangle K (in 3D the bubble is a quartic polynomial). Then the finite-dimensional subspaces for the pathogen and leukocyte , and nominal porosity Φ h ⊆ L 2 (Ω) are defined, respectively, as follows The pair (V h , Φ h ) is the well-known MINI-element, which is inf-sup stable in the context of saddle-point Stokes equations in their velocity-pressure formulation [28].
Then the fully discrete problem arises from (3.2) and for each time step n, we perform inner Newton-Raphson iterations from k = 0, . . . seeking (δc k+1 p,h , δc k+1 l,h , δu k+1 H h solutions to the unsymmetric matrix system where the entries δC k+1 p , δC k+1 l , δU k+1 , δP k+1 and δΦ k+1 f in the independent vector variable are the vectors containing all internal degrees of freedom associated with the discrete incremental solutions for all fields, and the operators in calligraphic letters appearing in the coefficient matrix and load vector from (5.2) are induced by the corresponding bilinear forms and linear functionals in (3.3). The functionals with superscript n on the right-hand side vector indicate that they also receive contributions from the backward Euler timediscretisation.

Schur complement based robust preconditioner
The numerical solution of problem (5.2) through direct methods is not feasible for large systems, which is a natural consequence of considering finer meshes for 3D geometries, needed to better capture anatomical details. This is the standard scenario in which Krylov space methods are the most useful, more specifically a GMRES method due to the indefinite and non-symmetric nature of the problem [69]. Still, scalable solvers require the construction of a robust preconditioner that allows for the Krylov iterations to remain roughly constant whenever (i) the amount of processors is increased and (ii) the mesh is refined. In this section we propose a preconditioner based on a two-level nested Schur complement inspired by the multiphysics nature of the model. If we consider a block matrix M given by the general structure with A invertible, a Gauss elimination procedure yields We note that if C is invertible, the same procedure can be applied with respect to it. Schur complement based preconditioners enjoy excellent theoretical properties, as the preconditioned system possesses at most 3 distinct eigenvalues [57], implying that it converges in at most 3 iterations of a Krylov subspace method, but this is seldom true in practice as the Schur complement block S = C − B 2 A −1 B 1 is computationally expensive to compute. Two complementary strategies to circumvent this problem are (i) to consider the full/lower triangular/upper triangular/diagonal part of (5.3) and (ii) to approximate S. In turn, two standard approaches for approximating S are the SIMPLE preconditioner given by S ≈ C − B 2 diag (A) −1 B 1 and the block diagonal approximation S ≈ C (for other possibilities and their comparison see, e.g., the comprehensive review [35]).
On the first level of the proposed preconditioner, we split the variables into poroelastic and chemotaxis, denoted by poro = (u, p, φ) and chem = (c p , c l ) respectively, and use this to write the linearised problem as where each block is given by

Using (5.3) we can write the inverse of (5.4) as
This idea has been applied in [33] for an FSI problem, and in [48] for the Boussinesq equations. The equations governing (c p , c l ) are of parabolic type and so their approximation does not pose a challenge, which indicates that the main difficulty lies in the poroelastic block approximation. Motivated by this, the proposed preconditioner is given by the following components: -On the first level we consider a lower-triangular Schur decomposition and on the second level we consider a full one. -The "chem" block is approximated by the action of the HYPRE-BoomerAMG [81] preconditioner. -Both blocks in the second level are approximated by an inexact GMRES solver preconditioned by an additive Schwarz method with a direct solver in each subdomain (using the MUMPS library [4]).
The inexactness is given by a relative tolerance of 0.1, and we highlight that an inexact solver in the subblocks gives rise to a preconditioner that changes from one iteration to the next one. This requires the use of a flexible GMRES algorithm (fGMRES) [68]. The use of an additive Schwartz preconditioner, or in general a domain decomposition one, alleviates the deteriorating scalability of AMG preconditioners for higher-order elements, in this case required by the displacement to satisfy an appropriate discrete inf-sup condition.
We close this section by noting that the choice of the type of Schur decomposition is not arbitrary. For this, we consider exact, lower triangular Schur complement preconditioners at both levels. In this way, the preconditioner adopts the following form After algebraic manipulations, it can be seen that the preconditioned system is given by where I .
The preconditioned system does not exhibit the expected property of having a small number of eigenvalues. This is mainly due to the use of a preconditioner instead of an exact inverse, where P poro = J poro . If such blocks coincide, then we would have which justified the use of an exact Schur preconditioner at the second level. This results in a preconditioned problem with only one eigenvalue.  Reference parameters of the model to be used in Section 6.1.

Computational experiments
We now turn to the presentation of numerical examples serving to illustrate the performance of the finite element scheme and to examine further the main features of the model. All routines have been implemented using the open source finite element library FEniCS [1]. A fixed tolerance of 10 −6 is used on the residuals for the convergence criterion of the Newton-Raphson iterative algorithm. We highlight that the error control of multiphysics problems is not trivial, for example in [16] as an error they use the norm of the increment, whereas in [80] they use the norm of the residual. We prefer the norm of the residual, as having an increment going to zero might mean stagnation and not convergence. Still, this can be improved by checking the error in each physics separatedly. A simplified sensitivity analysis varying one parameter at a time with respect to the reference parameters reported in Table 6.1 was performed. We limit the presentation of these results to only a few parameters and variables divided into three cases that highlight the couplings of the model. Young's modulus E, associated with the mechanical part, was investigated for the first case. In the second and third cases, associated with the immune response of the model, both phagocytosis rate, λ lp , and pathogen's reproduction rate, γ P , were studied. These parameters were chosen because they better highlighted the couplings in the model. They are also related to critical biological responses depending on their values. Changes in the value of γ P , for example, can be associated with pathogens with distinct proliferation abilities; changes in the value of λ lp can describe the ability of leukocytes in dealing with the invading pathogen; whereas changes in E represent tissues with different stiffness.
In all the three cases, the investigated parameter (E, λ lp , or γ P ) had its value doubled and halved with respect to the reference value reported in Table 6.1. To study these three scenarios, we assume that pathogens are responsible for triggering the inflammatory response. For this reason, the instant in which pathogens' concentration reaches its peak is used as a reference to collect data from the other variables of interest. Figure 6.1 shows the results of the sensitivity analysis. The time in which pathogens concentration reaches its maximum was distinct for each scenario only for the analysis of the influence of pathogen's reproduction rate (cf. Figure 6.1, third line), occurring after 12 days, six days, and 24 days, respectively, since the start The plots in the first row of Figure 6.1 presents the results for varying tissue stiffness. A stiffer tissue will imply a lower capacity to deform and accumulate fluid. With a smaller amount of fluid in the infected region, the concentration of pathogens and leukocytes in the tissue increases ( Figure 6.1, top left), i.e., the same amount of pathogens in the region is divided by a smaller fluid volume.
The second row of panels in Figure 6.1 presents the results for varying the phagocytosis rate, λ lp . A reduction in λ lp value can be related to a scenario in which the immune system is ineffective in dealing with the invading pathogen. This could be explained by immunodeficiency, for example. In this context, the concentration of pathogens will reach a higher and faster peak of infection for smaller values of λ lp . The lower capacity of phagocytosis by leukocytes leads to a considerable increase in the concentration of pathogens (Figure 6.1, centre left). The increase in pathogens' concentration induces an increase in the capillary permeability, which allows a more significant accumulation of fluid in the infected region ( Figure 6.1 centre middle). The accumulation of fluid in the region causes the pressure to increase. As a consequence of these events, the deformation suffered by the tissue is greater than the one observed for higher values of λ lp (Figure 6.1 centre right).
Note that an increase in γ P can be interpreted as a pathogen that can reproduce faster due to its intrinsic biological characteristics or due to the presence of favourable conditions (abundance of nutrients or temperature, for example). In fact, as the bottom row of Figure 6.1 shows, the peak of infection was higher and faster for larger γ P values, as expected. However, counter-intuitively, high values of γ P are associated with small deformations and small regions of oedema. As a result of more pathogens, the capillary permeability increases to allow more leukocytes to enter into the tissue (Figure 6.1, bottom left). The rapid increase in the concentration of leukocytes contained the infection quickly, preventing pathogens from remaining in the tissue long enough to spread. Since the time pathogens stay in the tissue is short, as well as the region in which they spread, the deformation suffered by the tissue also decreases ( Figure 6.1, bottom right).
In the third scenario, with the small value of γ P = 0.06, the peak of infection only occurs 24 days after the invasion begins, which leads to a weak inflammatory response. Thus, the concentration of leukocytes is reduced. The low leukocyte concentration results in a long time to remove pathogens that can spread over a larger region. The weak immune response results in more deformation when compared to the cases with faster pathogen dynamics, which highlights the couplings and strong nonlinearities in the model.
For this test we have also used the FEniCS environment along with the direct solver MUMPS.

Example 2: Compression and drainage of a poroelastic sample
With the aim of illustrating the coupling between poroelastic deformation and fluid flow using the formulation (2.2a)-(2.2c) and in the absence of the immune system chemotaxis, we simulate the incremental compression of a box Ω = (0, 1) 2 m 2 . The bottom edge of the boundary constitutes Γ (where according to (2.6a), the deformation is zero and we impose zero fluid flux), whereas Σ is conformed by the top edge (where a sinusoidal-in-time distributed traction t Σ = 2000 sin(πt)n is applied and a constant fluid pressure p in = 0.2 MPa is considered) as well as the vertical walls (where we prescribe zero traction and zero fluid pressure). We employ a non-homogeneous initial Lagrangian porosity that is not symmetric with respect to the centre of mass of the domain, given by φ 0 (x) = 3 5 + 1 10 sin(x 1 x 2 ), the isotropic power-law porosity-dependent permeability, a compressible neo-Hookean strain energy density yielding the nominal effective stress S eff = µ s (F t − F −1 ) + λ s lnJF −1 , and the parameters The scheme characterised by the finite element spaces in (5.1) is used on a crossed-shaped uniform mesh with 10000 cells, and the coupled system is simulated until t final = 1 s using a constant time step ∆t = 0.1 s. The system is considered initially at rest (u(0) = 0, p(0) = 0, φ f (0) = φ 0 ). This test also serves to assess the robustness of the method. We run a set of simulations varying the Poisson ratio from 0.2 to 0.499999. Regardless of the parameter regime, the numerical results do not show spurious oscillations of pressure, nor unexpectedly small displacements or nonphysical distortions or checker-board patterns in the porosity field. In Figure 6.2 we plot deformed geometries and field variables showing the progressive compression of the porous block over three time instants, and using ν = 0.33 and ν = 0.495. As expected, the fluid pore pressure does not have a uniform distribution, and it moves toward Γ on the bottom edge.
In addition, we perform a benchmark test where again a time-harmonic loading defined by t Σ = 1 5 (λ s + 2µ s ) sin( 2π 15 t)n is applied, now only on part of the top edge of the box Ω = (0, 8) × (0, 5) m 2 (on the segment 0 ≤ x 1 ≤ 1 m, x 2 = 5 m). The setup of this validation example proceeds similarly as in [49,Sect. 4.2] (see also [84,Sect. 4.2] and [51,65]). The boundary conditions in these references permit drainage (fixing and we now use the normalised and isotropic Kozeny-Carman porosity-dependent permeability. For this test the mesh is unstructured and graded near the top-left corner. Figure 6.3 shows snapshots of the approximated field variables at times t = 4, 12, 16 s and the bottom panels also portray transients of the vertical displacement and porosity recorded at two locations (node A≈(0.5,4.5), directly below the centre of the footing strip and node B≈(8.0,4.5) approximately at the same height but on the right wall). These plots show the same qualitative behaviour encountered in, e.g., [51].

Example 3: Errors with respect to manufactured solutions
Since there is no analytical solution currently available for the coupled system (2.2a)-(2.4b), the accuracy of the finite element discretisation is investigated using the following closed-form manufactured solutions conditions are considered for the displacement and fluid pressure. On the left side of the boundary Γ , we prescribe the exact displacement from (6.1) and an exact traction Pn on the remainder of the boundary, Σ. The fluid pressure is constrained to match the exact one on Σ and an exact fluid pressure flux is imposed on Γ .
We use the same neo-Hookean material law as in Example 1, and employ the isotropic Kozeny-Carman permeability κ KC iso . The following parameter values are selected they are all regarded adimensional and do not have physical relevance in this case, as we will be simply testing the convergence of the finite element solutions.
We generate successively refined simplicial grids and use a sufficiently small (non dimensional) time step ∆t = 0.01 and simulate a sufficiently short time horizon t final = 3∆t, to guarantee that the error produced by the time discretisation does not dominate. Errors between the approximate and exact solutions are plotted against the number of degrees of freedom in Figure 6.4. This error history plot confirms the optimal convergence of the finite element scheme (in this case, first order) given by the Remark 4.1, for all variables in their respective norms, where a slightly better rate, of about 1.3, is seen for the porosity. In addition, the average iteration count for the Newton-Raphson method (for each refinement level and for all time steps) is five.

Example 4: Localisation of oedema regions
In the next test we present and discuss the coupling dynamics of the model (2.2)-(2.4), considering a twodimensional domain Ω = (0, 4) × (0, 4) cm 2 , and a time step size of ∆t = 0.1 days. The parameters described in the Table 6.1 were used. Boundary conditions were applied as follows: u = 0 on the left edge, and p = 0 for the right, top and bottom edges of the domain, and the following initial conditions were used: c p,0 = 0.001 in the region (x 1 − 2) 2 + (x 2 − 2) 2 ≤ 0.03 and c p,0 = 0.0, otherwise. Also, we consider c l,0 = 0.003, φ 0 = 0.2 and p 0 = 0 in the entire domain. Figure 6.5 presents the results of a localised oedema formation. Each row presents snapshots of the variables c p , c l , p, u and φ f , respectively, at three selected time instants.
The dynamics of the immune response can be observed in Figure 6.5 in terms of c p and c l . After the pathogens appear in the tissue, they start to grow (Figure 6.5 at time t = 8 days) and reach their peak concentration ( Figure 6.5 at time t = 13 days). However, in response to these events, leukocytes migrate to the infected site and their concentration also grows ( Figure 6.5, c l column). As a consequence of the presence of leukocytes, pathogen concentration starts to decrease, as Figure 6.5 shows for c l at time t = 18 days.
The presence of pathogens triggers a mechanical response due to a sequence of couplings in the model. The endothelium increases its permeability to allow leukocytes to leave the bloodstream and enter the tissue. This increase in the endothelium permeability in turn increases interstitial fluid and pressure ( Figure 6.5 for φ f and p at time t = 8 days). When the concentration of pathogens reaches its peak ( Figure 6.5 for c p at time t = 13 days), it is possible to observe that the amount of fluid phase φ f in the tissue also reaches its maximum value. An increase of liquid in the region leads to an increase in both pressure and displacement fields ( Figure 6.5 for u and p at time t = 13 days). Finally, when pathogens c p are almost eliminated by the leukocytes at time t = 18 days, the dynamics of pressure, fluid phase fraction and displacements, tend to return to their initial values. 6.5 Example 5: Immune system dynamics and poro-hyperelasticity in a left ventricle Next we conduct a series of tests using a ventricular geometry segmented from patient-specific images [78], and where synthetic muscle fibre and collagen sheet directions are constructed using a Laplace-Dirichlet  rule-based method [66]. On the basal surface we prescribe zero normal displacement and zero fluid pressure flux, the epicardium is considered stress-free and with a prescribed fluid pressure, and on the endocardium we impose a time-dependent traction t Σ = m 0 sin 2 (πt)n with m 0 = 0.1 N/cm 2 , and an endocardial fluid pressure p endo = p 0 sin 2 (πt) having the same time period. No-flux conditions are used for the concentrations on the whole boundary.
For this example we use the Holzapfel-Ogden constitutive strain energy stated in (2.1), and employ the anisotropic Kozeny-Carman permeability κ KC introduced in (2.3); while the remaining model parameters used for the 3D ventricular test assume the values λ s = 27.293 kPa, µ s = 3.103 kPa, D p = 3 × 10 −3 /φ 0 cm 2 /h, π i = 10 mmHg, n = 5, α = 0.5, An initial concentration of pathogens is considered on a transmural strip, and we simulate the dynamical behaviour of the coupled system over several minutes. The evolution of the pathogens and leukocytes distribution on the deformed ventricular geometry is shown in Figure 6.6, where we also show fibre distribution and mechanical fields at half-time. One can observe that the presence of pathogens induces the local entry of leukocytes. Fluid also enters the interstitial space of the tissue, increasing pressure. Although leukocytes can remove part of the pathogens, some of them diffuse through almost the entire domain. As this pathogens' Fig. 6.6: Example 5. Assigned fibre distribution streamlines, displacement magnitude, fluid pressure distribution, and porosity on a patient-specific left ventricular geometry at t = 14 (top), and evolution of pathogens and leukocytes concentration (middle and bottom rows, respectively, seen from a slightly different angle).
wave sweeps a significant part of the left ventricle, more parts of the heart are impacted by the inflammatory response and, as a consequence, diffuse oedema is formed. Numerical simulations show that the propagation of the front of pathogens likely depends on the diffusion and replication of pathogens and the wave tail likely depends on the diffusion and efficiency of the leukocytes [52]. The numerical results indicate that, although fluid phase and pressure increase due to the presence of the pathogen, they result in small changes in displacement. In fact, new diagnostic tools have been studied to increase the non-invasive diagnostic accuracy of myocarditis and other myocardial pathologies using for this purpose the detection of increased water contents at some regions instead of an increase in tissue displacement. See, for example, [47,72,73].
To conclude this section, for the ventricular geometry we proceed to briefly test the preconditioner described in Section 5.2. To measure its performance, we look at the total and average number of Krylov iterations in two different meshes corresponding to the left ventricle geometry and with respect to P 1 and P b 1 (bubble enriched element) conforming finite elements for the displacement to study the impact of the inf-sup Table 6.2: Example 5. Iteration counts, total and (average), for fGMRES with the proposed preconditioner, and for either continuous and piecewise linear approximations of displacement, or piecewise linear displacements with bubble enrichment. Note that an increase in the total iterations, for constant average iterations, occurs due to a different number of Newton iterations. stability in the performance. The results are shown in Table 6.2, and as we have observed that numerically the most difficult time instant is the first one, we report the performance only on it.
We highlight that the number of average Krylov iterations is roughly constant with an increasing number of cores. This can be expected due to the use of inexact solvers in the sub-blocks of the Schur complement preconditioner, which yield an adequate approximation of the exact Schur preconditioner, which has at most 3 distinct eigenvalues. It is particularly interesting to observe that this performance is obtained by using only the action of an AMG preconditioner for the chemotaxis block, meaning that (i) the coupling between the chemotaxis and poroelastic models is not very strong in passive cardiac simulations and that (ii) such block is not necessarily computationally challenging.

Final remarks
We have studied a general model capturing the phenomenological features of the interaction between chemotaxis of the immune system in saturated poroelastic media admitting large deformations. The essential properties of the new model include a formulation in terms of displacement, fluid pressure, and Lagrangian porosity, coupled with concentrations for leukocytes and pathogens. In particular, the poro-hyperelastic compartment of the model can be identified, after linearisation and adequately choosing the initial guess, with the three-field formulation for Biot's poroelasticity from [60,77]. We have proposed a finite element method together with splitting strategies, and have studied their properties in terms of accuracy and iteration count. The realisation of the coupling is general enough to accommodate different types of model specifications, including diverse hyperelastic solid laws and showing the interplay of solid deformation, effective stress, and pore fluid pressure build-up. The model also includes species transport through the total velocity, and the system is solved by means of a monolithic finite element method in saddle-point form.
Some extensions in the development of this work include the study of long-term behaviour of the oedema formation, as well as performing a thorough exploration of different effects from introducing lymphatic sinks for c p and c l . We also plan to incorporate active stress due to calcium release (in turn generated by the immune reaction-diffusion) and assess quantitative differences in the classical diffusion case. In a more applicative context, a more comprehensive sensitivity analysis is still required to determine the most relevant coupling mechanisms. Also, the validation and verification of the model against patient-specific data are still to be conducted. We have expanded the application discussed in Section 6.5 to address further modelling and numerically oriented investigations in the recent work [52]. Regarding numerical schemes, a possible next step is to use mixed formulations for the immune system equations (following, e.g., [43]) to produce mass conservative discretisations.
Further investigation is necessary, for instance, regarding the specific form of the anisotropic porosity as well as in designing new coupling mechanisms that will contribute to a better understanding of the formation and termination of myocarditis and myocardial oedema.