### Patient data and image processing

Cardiac MRI was performed in a 27-year-old male patient (85*kg*, 171*cm* height) after 5 days of the beginning of symptoms. The patient had undergone a previous invasive angiography that revealed no coronary artery disease. MRI was performed on a 3.0T scanner (Siemens Skyra, Syngo version E11C, Erlangen, Germany) with a 30-channel phased array coil. T2 parametric maps were obtained in a short and 4-chamber long axis using a free breathing sequence with T2-prepared single-shot steady-state free precession (SSFP) imaging [52]. LGE images were obtained after the injection of 0.2*mmol* of gadolinium-based contrast using a sequence with these parameters: gradient-echo inversion recovery (GRE-IR) sequence (TR 6.0*ms*; TE 4.18*ms*; matrix 256×128; FOV 350×260; resolution 1.37×2.0×8*mm*; FA 25o; BW 130*Hz*/*px*; TI 400*ms*; 25 segments; 2 RR intervals; 9 images were obtained in short axis, 3 images in two-chamber view and 3 images in four-chamber views [53]). Figure 1a presents the LGE results.

Post processing was performed using a dedicated workstation with specialized software (CVI42 v5.9, Circle CVI, Calgary, Canada). All images were analyzed by the same reader (JLF, ten years of experience in CMR). For T2 mapping a single region-of-interest was drawn in the septum and lateral wall to obtain the T2 data (Fig. 1b). For LGE images, visual qualitative analysis was performed, classifying the findings as positive or negative for LGE without automated quantification.

### Mathematical model

The mathematical model used to describe the pathophysiology of edema formation was proposed in a previous work [49]. This model includes the inflammatory part which describes the interaction between a pathogen and the human immune system. Along with the inflammatory dynamics, this model also uses a poroelastic model, based on Biot poroelasticity theory, to describe the mechanical deformation of the simulated tissue. Afterwards, this model was proposed to describe the effects of infectious myocarditis in the heart [50].

The model couples the immune response to an extracellular edema formation. The model considers a fluid-saturated porous medium *Ω* to represent the tissue and the pathogen is modelled by:

$$ \left\{\begin{array}{ll} \frac{\partial (\phi_{f} C_{p})}{\partial t} = \nabla \cdot (D_{b} \nabla C_{p}) -r_{b} + q_{b} & \ \text{ in}\ \Omega\times I,\\ D_{b}\nabla C_{p} \cdot \boldsymbol{n} = 0 & \ \text{ in}\ \partial \Omega \times I,\\ C_{p}\left(\boldsymbol{x},0\right)=C_{p0}(\boldsymbol{x}) & \ \text{ in}\ \Omega, \end{array}\right. $$

(1)

where \(\Omega \subset \mathbb {R}^{2}\) and \(I=(0,t_{f}]\subset \mathbb {R}^{+}\) is the time interval, *n* is the normal vector outward the boundary domain *∂**Ω*, \(C_{p}:\Omega \times I\rightarrow \mathbb {R}^{+}\) is the pathogen concentration in interstitial fluid, *ϕ*_{f} is the porosity, *D*_{b} is the diffusion coefficient of the pathogen through the interstitial fluid, *q*_{b} denotes the pathogen reproduction, and *r*_{b} denotes the pathogen death due to leukocytes action. Diffusion is defined as the spread of particles from regions of higher concentration to regions of lower concentration. The terms *q*_{b} and *r*_{b} represent the source of pathogens and phagocytosis of pathogens by leukocytes, respectively, which are defined by:

$$\begin{array}{*{20}l} q_{b} &= c_{p} C_{p}, \end{array} $$

(2)

$$\begin{array}{*{20}l} r_{b} &= \lambda_{nb} C_{l} C_{p}, \end{array} $$

(3)

where *c*_{p} is the pathogen growth rate in the interstitial tissue; *C*_{l} is the leukocytes concentration in the interstitial fluid and *λ*_{nb} is the leukocytes phagocytosis rate [51].

The leukocyte differential model is represented by:

$$ {{}\begin{aligned} \left\{\begin{array}{ll} \frac{\partial (\phi_{f} C_{l})}{\partial t}= \nabla \cdot (D_{n} \nabla C_{l} - \chi_{nb} C_{l} \nabla C_{p}) - r_{n} + q_{n} &\text{ in}\ \Omega\times I,\\ (D_{n}\nabla C_{l} -\chi_{nb} C_{l} \nabla C_{p})\cdot \boldsymbol{n} = 0 & \text{ in}\ \partial \Omega \times I,\\ C_{l}\left(\boldsymbol{x},0\right)=C_{n0}(\boldsymbol{x}) & \text{ in}\ \Omega, \end{array}\right. \end{aligned}} $$

(4)

where \(C_{l}:\Omega \times I\rightarrow \mathbb {R}^{+}\) is the leukocytes concentration in the interstitial fluid, *D*_{n} is the diffusion coefficient of the leukocyte through the interstitial tissue, *χ*_{nb} is the leukocyte chemotaxis rate, *q*_{n} is the leukocyte source, which represents the leukocyte extravasation from bloodstream to interstitium, and *r*_{n} denotes the leukocyte death due to both apoptosis and pathogen phagocytosis. These terms are defined by:

$$\begin{array}{*{20}l} q_{n} &= \gamma_{n} C_{p} (C_{n,max} - C_{l}), \end{array} $$

(5)

$$\begin{array}{*{20}l} r_{n} &= \lambda_{bn} C_{l} C_{p} + \mu_{n}C_{l}, \end{array} $$

(6)

where *C*_{n, max} is the leukocytes concentration in the bloodstream, and *γ*_{n} is the leukocyte permeability to capillaries microvascular wall [51]. Here, *λ*_{bn} represents the death rate after phagocytosis and *μ*_{n} represents the natural leukocyte decay, once some of them are preprogrammed to die a short time after they leave the bloodstream [8].

With respect to the mechanical modeling, the stresses related to the tissue can be described by two parts: one which is caused by the hydrostatic pressure of the interstitial fluid filling the pores, and the other caused by the average stress in the solid part, i.e. the extracellular matrix and tissue cells.

Once *Ω* is a fluid-saturated porous medium representing the tissue, therefore a 2-phase system of equations derived from the principle of mass conservation is considered. In this system of equations, one of the equations can be used to describe the fluid (f) and the other the solid (s) phase. The fluid phase is mainly composed of blood plasma, whereas the solid phase represents the extracellular matrix along with the tissue cells. This system is given by:

$$ \left\{\begin{array}{ll} \frac{\partial (\phi_{f} \rho_{f})}{\partial t} + \nabla \cdot \left(\phi_{f} \rho_{f} \boldsymbol{v_{f}} \right) = q & \text{for the fluid phase,}\\ \frac{\partial (\phi_{s} \rho_{s})}{\partial t} + \nabla \cdot \left(\phi_{s} \rho_{s} \boldsymbol{v_{s}} \right) = 0 &\text{for the solid phase,} \end{array}\right. $$

(7)

where *q*_{c} and *q*_{l} represent the capillary bed network and the lymphatic system influence in *Ω*, respectively, and *q*=(*q*_{c}+*q*_{l})*ρ*_{f}. In addition, *ϕ*_{α}, *ρ*_{α} and *v*_{α} are porosity, density and velocity for each phase, where *α*=*f* for fluid and *α*=*s* for solid.

The plasma flux of the lymphatic capillaries is influenced by the interstitial fluid pressure [54]. The influence of the lymphatic system on the interstitial fluid dynamics, described by the term *q*_{l} in Eq. (7), is given by:

$$ q_{l}(P)= -\bigg(q_{0}\bigg(1+ \frac{V_{max} (P-P_{0})^{n}}{K_{m}^{n} + (P-P_{0})^{n}}\bigg)\bigg), $$

(8)

where *q*_{0} is the normal lymph flow, *P*_{0} is the normal interstitial fluid pressure. Moreover, *V*_{max} is the maximum lymph flow, *K*_{m} is the half-live, i.e. the value of *P* which the lymph flow corresponds to \(\frac {V_{max}}{2}\), and *n* is the Hill coefficient. These parameters are normally obtained by experimental data [55].

Following our previous work [43], the capillary bed network term *q*_{c} is given by:

$$ q_{c}(P) = c_{f}(P_{c} - P -\sigma(\pi_{c} - \pi_{i})), $$

(9)

where *c*_{f} is the filtering coefficient given by *L*_{p}(*S*/*V*), *L*_{p} is the hydraulic permeability of the microvascular bed wall, and (*S*/*V*) is the vessel surface area per unit volume; *P* and *P*_{c} are the fluid pressure in the interstitium and the capillary pressure, respectively; *π*_{c} and *π*_{i} are the capillary oncotic pressure and the interstitial oncotic pressure due to the plasma protein, respectively; and the reflection coefficient to plasma proteins is denoted by *σ*∈[0,1]. Eq. (9) is known as Starling equation [56] and is widely used [42, 54, 57].

In order to couple the influence of the inflammatory reaction to the fluid phase, the following model [43] for the hydraulic permeability and the oncotic reflection coefficient was used:

$$ L_{p}(C_{p}) = L_{p0}(1+c_{bp}C_{p}), $$

(10)

where *L*_{p0} is the health tissue hydraulic permeability of the capillary wall, *C*_{p} is the pathogen concentration, and *c*_{bp} is the influence of the pathogen infection in the hydraulic permeability.

The coupling of oncotic reflection coefficient with the inflammatory response is given by the following relation:

$$ \sigma(C_{p}) = \frac{\sigma_{0}}{(1+c_{br}C_{p})}, $$

(11)

where *σ*_{0} is the oncotic reflection coefficient in a non-inflamed tissue, *C*_{p} is the pathogen concentration, and *c*_{br} is the influence of pathogens in the reflection coefficient.

Assuming \(\boldsymbol {v_{s}} = \frac {\partial \boldsymbol {U}}{\partial t}\), then Eq. (7) can be rewritten as follows:

$$ \left\{\begin{array}{ll} \frac{\partial (\phi_{f} \rho_{f})}{\partial t} + \nabla \cdot \left(\phi_{f} \rho_{f} \boldsymbol{v_{f}} \right) = q & \text{for the fluid phase,}\\ \frac{\partial (\phi_{s} \rho_{s})}{\partial t} + \nabla \cdot \left(\phi_{s} \rho_{s} \frac{\partial \boldsymbol{U}}{\partial t} \right) = 0 &\text{for the solid phase.} \end{array}\right. $$

(12)

Considering that both fluid and solid densities are constant and that *ϕ*_{s}+*ϕ*_{f}=1, then Eq. (12) can be rewritten to obtain the following system of equations:

$$ \left\{\begin{array}{l} \nabla \cdot \boldsymbol{v_{D}} + \nabla \cdot \frac{\partial \boldsymbol{U}}{\partial t} = q_{c} + q_{l}, \\ \boldsymbol{v_{D}} = -\frac{\boldsymbol{K}}{\mu}\nabla P, \end{array}\right. $$

(13)

where *v*_{D} is the velocity of the fluid phase under a porous media approach, *K* is the porous media permeability, and *μ* is the fluid viscosity. The second equation is based on Darcy’s Law, and is often referred to as Darcy velocity [58].

The closure of the PDE system given by Eq. (13) comes when the momentum conservation is considered for each phase, resulting:

$$ \left\{\begin{array}{ll} \nabla\cdot(\boldsymbol{\sigma_{s}}) + \boldsymbol{\hat{T_{s}}}= \boldsymbol{0} & \text{for the solid phase,}\\ \nabla\cdot(\boldsymbol{\sigma_{f}}) + \boldsymbol{\hat{T_{f}}}= \boldsymbol{0} & \text{for the fluid phase,} \end{array}\right. $$

(14)

where *σ*_{f} and *σ*_{s} are the fluid and solid phases tensors, respectively, and \(\boldsymbol {\hat {T_{f}}}\) and \(\boldsymbol {\hat {T_{s}}}\) are the interaction forces between the phases.

Summing up the equations of each phase and considering \(\boldsymbol {\hat {T_{f}}} + \boldsymbol {\hat {T_{s}}} = 0\) results in the following equation for the mixture:

$$ \nabla \cdot (\boldsymbol{\sigma_{f}} + \boldsymbol{\sigma_{s}}) = \boldsymbol{0}. $$

(15)

The solid phase is considered as an isotropic elastic porous medium, i.e. a Hookean material model [45, 58]. Therefore, the solid phase stress tensor is given by:

$$ \boldsymbol{\sigma_{s}} = \lambda_{s} (\nabla \cdot \boldsymbol{U}) \boldsymbol{I} + 2\mu_{s} \boldsymbol{\varepsilon}(\boldsymbol{U}), $$

(16)

where \(\varepsilon (\boldsymbol {U}) = \nabla \boldsymbol {U}^{s} = \frac {1}{2}(\nabla \boldsymbol {U} + \nabla \boldsymbol {U}^{T})\), *I* is the identity matrix, and *λ*_{s} and *μ*_{s} are the Lamé parameters. The fluid phase tensor can be described according to the following relation [45, 58]:

$$ \boldsymbol{\sigma_{f}} = - P \boldsymbol{I}, $$

(17)

where *P* is the hydrostatic pressure at the interstitium and *I* is the identity matrix.

Finally, rearranging Eqs. (13) and (15), the following system of equations can be obtained in terms of *P* and *U*:

$$ \left\{\begin{array}{l} (\lambda_{s} + \mu_{s})\nabla (\nabla \cdot \boldsymbol{U}) + \mu_{s} \nabla^{2} \boldsymbol{U} - \nabla P = \mathbf{0},\\ \nabla \cdot \left(\boldsymbol{k} \nabla P\right) = \frac{3}{3\lambda_{s} + 2\mu_{s}} \frac{\partial P}{\partial t} + q, \end{array}\right. $$

(18)

where \(\boldsymbol {k} = \frac {\boldsymbol {K}}{\mu }\) is the mobility tensor and *U* is the displacement field.

The influence of pressure \(P:\Omega \rightarrow \mathbb {R}\) is governed by the following problem:

$$ \left\{\begin{array}{ll} \nabla \cdot \left(\boldsymbol{k} \nabla P\right) = \frac{3}{3\lambda_{s} + 2\mu_{s}} \frac{\partial P}{\partial t} + q & \text{ in}\ \Omega \times I, \\ \boldsymbol{k} \nabla P \cdot \boldsymbol{n} = 0 & \text{ on}\ \partial \Omega \times I,\\ P\left(\boldsymbol{x},0\right)=P_{0}(\boldsymbol{x}) & \text{ in}\ \Omega, \end{array}\right. $$

(19)

whereas the governing equation for the displacement field is given by:

$$ \left\{\begin{array}{ll} (\lambda_{s} + \mu_{s})\nabla (\nabla \cdot \boldsymbol{U}) + \mu_{s} \nabla^{2} \boldsymbol{U} - \nabla P = \mathbf{0} & \text{ in}\ \Omega,\\ \boldsymbol{U} = \boldsymbol{0} & \text{ on}\ \partial \Omega,\\ \end{array}\right. $$

(20)

where \(\boldsymbol {U} : \Omega \rightarrow \mathbb {R}^{2}\); *∂**Ω* defines the domain boundary which, in this case, it was applied a homogeneous Dirichlet boundary condition.

### Computational model

To obtain a numerical approximation for all the PDEs presented before we used the finite element method (FEM) due to its flexibility in handling complex geometries and boundaries such as long axis slice of a heart. Here, the FEniCS software library [59, 60] was used to obtain a numerical approximation of the variational formulation of the problems considered in this work. In particular, the models were discretized using linear triangles for spatial discretization and time-derivatives were approximated by a backward difference. The leukocyte model was solved using the SUPG (streamline upwind Petrov-Galerkin) form of the finite element method to ensure stability due to the presence of the convective term [61].

### Mesh generation

The software Gmsh [62] was used to segment the cardiac MRI exam and generate the finite element mesh which was used to perform the simulations with FEniCS. The geometry description was generated based on Fig. 1b allowing to create a patient-specific 2D finite element mesh with 5015 triangles. Figure 5a shows the resulting mesh. In addition, Gmsh was used to generate a binary image (Fig. 4a), based on the edema area identified by the specialist (Fig. 1b). The idea is to make easier the comparison between the patient edema area and the simulated one (presented in Fig. 4b).

### Simulations setup

Tables 1 and 2 present all parameters values used in Eqs. (1) and (4), respectively. These values were based on previous studies [43, 44, 49, 63]. The values used in Eqs. (19) and (20) are shown in Table 3. To perform the simulations, it was also necessary to set up proper boundary and initial condition of Eqs. (1), (4), (19) and (20). All their functions and values are shown in Table 4.

The patient-specific mesh was created after the segmentation of the long axis cardiac MRI image of the patient (see Methods and Fig. 5a). Lymph vessels represent about 2.9*%* of tissue [64]. To perform the simulations, the finite element mesh was randomly and uniformly filled with lymph vessels, of the size of one element each, with a volume fraction of 2.9*%*. The remaining elements are considered to be only under blood capillary influence. Figure 5b presents an example of a distribution of lymph vessels in the domain, where triangular elements representing lymph vessels are marked in white.