Skip to main content

A personalized computational model of edema formation in myocarditis based on long-axis biventricular MRI images



Myocarditis is defined as the inflammation of the myocardium, i.e. the cardiac muscle. Among the reasons that lead to this disease, we may include infections caused by a virus, bacteria, protozoa, fungus, and others. One of the signs of the inflammation is the formation of edema, which may be a consequence of the interaction between interstitial fluid dynamics and immune response. This complex physiological process was mathematically modeled using a nonlinear system of partial differential equations (PDE) based on porous media approach. By combing a model based on Biot’s poroelasticity theory with a model for the immune response we developed a new hydro-mechanical model for inflammatory edema. To verify this new computational model, T2 parametric mapping obtained by Magnetic Resonance (MR) imaging was used to identify the region of edema in a patient diagnosed with unspecific myocarditis.


A patient-specific geometrical model was created using MRI images from the patient with myocarditis. With this model, edema formation was simulated using the proposed hydro-mechanical mathematical model in a two-dimensional domain. The computer simulations allowed us to correlate spatiotemporal dynamics of representative cells of the immune systems, such as leucocytes and the pathogen, with fluid accumulation and cardiac tissue deformation.


This study demonstrates that the proposed mathematical model is a very promising tool to better understand edema formation in myocarditis. Simulations obtained from a patient-specific model reproduced important aspects related to the formation of cardiac edema, its area, position, and shape, and how these features are related to immune response.


According to the World Health Organization (WHO) and International Society and Federation of Cardiology (ISFC), myocarditis is defined as an inflammatory disease of the heart muscle [1, 2]. Recent studies address myocarditis in postmortem diagnosis of sudden cardiac death ranging from 2% to 44% in young adults [35]. Another evidence shows that subclinical myocarditis can be the cause of ventricular fibrillation: 42% of sudden cardiac deaths that occurred among US army recruits were associated with myocarditis [3]. The causes of myocarditis can be divided into infectious, immune-mediated and toxic. This paper focuses on infectious myocarditis which can be bacterial, spirochaetal, fungal, protozoal, parasitic rickettsial and viral [5].

When a living body is affected by an injury or is infected by an antigen, a complex chain of events is triggered, known as inflammatory response. Inflammation is characterized by edema, redness, warmth and pain. All of them emerges from complex physiological processes such as vasodilatation, increase of capillary permeability and pathogen phagocytosis [68].

Inflammatory edema is one consequence of an immunological response to an invading pathogen. The immune cells reach the infected area due to an increased vascular permeability, which also allows an increased plasma influx, and fluid accumulation. This process is known as extravasation and is an important mechanism to the edema formation [6].

Edema may be found in several tissues of the body such as lungs [9, 10], hands, arms, legs [1113] and heart [4, 14, 15]. Heart edema is of great research interest because it is commonly observed in myocarditis and is one of the criteria used in the medical practice to its diagnosis [4, 16]. Although it is invasive and has significant sampling error, endomyocardial biopsy (EMB) is considered the “gold standard” method for detecting myocardial inflammation and interstitial fibrosis. New diagnostic tools have been studied to increase non-invasive diagnostic accuracy of myocarditis and other myocardial pathologies. Cardiovascular magnetic resonance (CMR) imaging provides non-invasive assessment of cardiac function, anatomy and tissue characterization. New techniques have been developed and T2 parametric mapping allows the detection of myocardial edema due to increased water content at the edema regions with greater accuracy [1719].

Mathematical modeling is a extremely useful tool to uncovers complex physiological dynamics and interactions of biomedical phenomena and systems [2023]. In particular, distinct approaches have been used to model the interaction between the human immune system (HIS) and a pathogen, such as differential [21, 24, 25] and stochastic [26, 27] equations, cellular automata [28, 29], and agent-based models (ABM) [3032]. Models based on differential equations are widely used is the field. Ordinary Differential Equations (ODEs) are adopted when only the evolution of HIS populations over time is important, whereas Partial Differential Equations (PDEs) are used to capture spatiotemporal dynamics. Several studies have used ODEs [21, 24, 33, 34] and PDEs [25, 3539] to model the HIS.

To capture edema formation, the mechanical behavior of tissue and fluids also needs to be considered. Models for Interstitial fluid pressure (IFP) dynamics can be found in the literature [4042]. One important model that combines IFP dynamics and immune response was first suggested in our previous studies [43, 44]. To include tissue deformation, the classical theory of poroelasticity mechanics of a fluid-saturated porous media that was first proposed in 1941 [45] has been widely used in several studies [4648]. Two recent works proposed general models in one spatial dimension [49] and a two-dimensions [50] that couple the immune system response with poroelasticity theory. The models were used to describe inflammatory edema formation by also including the modeling of fluid accumulation.

In this work, we extend previous models of edema formation [49, 50] by developing and verifying a patient-specific model for edema formation in myocarditis in a two-dimensional domain. For this purpose, cardiac MRI images were obtained from a patient with myocarditis and were used to develop a specific geometric model of the patient’s heart. In addition, T2 parametric mapping was calculated from MRI images obtained from the same patient to quantitatively describe the edema in terms of location, shape and area. Patient-specific simulations using the new proposed hydro-mechanical mathematical model allowed us to correlate spatiotemporal dynamics of representative cells of the immune system, such as leucocytes and the pathogen, with fluid accumulation and cardiac tissue deformation. Finally, important features of the simulated edema, such as area, location and shape quantitatively matched the clinical ones extracted from the T2 mapping.


Patient data

Figure 1 presents Late Gadolinium Enhanced (LGE) image and T2 parametric map obtained in a 4-chamber long axis. LGE image shows the presence of an area of enhancement in the mid-wall portion of the lateral wall of the left ventricle (arrowheads). The localization of these areas of LGE are typical for patients with myocarditis and represent one of the main criteria for its diagnostic by MRI. T2 parametric map at the same 4-chamber slice position quantifies the edema. The region-of-interest 1 (ROI-1) was created using the location indicated by the LGE image. The region-of-interest 2 (ROI-2) was chosen far from this location. In the ROI-1 a T2 value of 44.8ms quantitatively identifies an edema when compared to the T2 value of 34.5ms in ROI-2. Therefore, the increased T2 value identified in the lateral wall colocalized with the LGE area seen in Fig. 1a suggests the presence of edema in comparison to an apparently normal area in the septum (ROI-2). The area of the edema measured by cardiac MRI is equal to 1.755cm2.

Fig. 1
figure 1

a LGE 4-chamber images showing the presence of an area of enhancement in the mid-wall portion of the lateral wall of the left ventricle (arrowheads). Heterogeneous, focal areas of LGE can also be identified in the septal wall. The localization of these areas of LGE are typical for patients with myocarditis and represent one of the main criteria for its diagnostic by MRI. b T2 parametric maps of the same patient in A, at the same 4-chamber slice position. The region-of-interest 1 in the lateral wall shows a T2 value of 44.8ms versus a T2 in the septum (ROI 2) of 34.5ms. The increased T2 identified in the lateral wall colocalizes with the LGE area seen in A and suggests the presence of edema/inflammation in comparison to an apparently normal area in the septum

Numerical results

Figures 2 and 3 show the simulation result considering Eqs. (1), (4), (19) and (20) after the simulation of 5 and 10 h of infection, respectively. The results were obtained solving the system of PDEs presented in section Methods using as parameters, initial and boundary conditions the values presented in the previous section. Four results are presented: pathogen (panel A) and leukocytes concentration (panel B), pressure (panel C) and displacements (panel D). Although we consider an isotropic elastic tensor, the displacements are larger in the y-axis than in the x-axis (Fig. 3d). This occurs due to the homogeneous Dirichlet boundary conditions imposed on the boundary of the domain (both on the epicardium and on the endocardium). In addition, the magnitudes of the displacements become greater in the border of the infection than in its center. This happens since the points in the vicinity of the infected site are moving in relation to the central points.

Fig. 2
figure 2

Simulation results of the edema formation using Eqs. (1), (4), (19) and (20) considering the model parameters presented in Tables 1, 2 and 3 with initial and boundary conditions defined in Table 4 after 5h of simulation. Panels A and B show the interactions between a pathogen infection and leukocytes considering an initial infection on edematous epicardium, while panels C and D show its consequences to interstitial fluid pressure and the displacement field

Fig. 3
figure 3

Simulation results of the edema formation using Eqs. (1), (4), (19) and (20) considering the model parameters presented in Tables 1, 2 and 3 with initial and boundary conditions defined in Table 4 after 10h of simulation. Panels a and b show the interactions between a pathogen infection and leukocytes considering an initial infection on edematous epicardium, while panels c and d show its consequences to interstitial fluid pressure and the displacement field

Table 1 Parameters values for Eq. (1) based on [43, 44, 49, 63]
Table 2 Parameter values for Eq. (4) based on [43, 44, 49, 63]
Table 3 Parameter values for Eqs. (20) and (19) based on [43, 44, 49, 63]
Table 4 Initial and boundary condition

Figure 4a shows the patient edema and Fig. 4b shows the numerical result of the fluid phase obtained by Eq. (7). An isoline is placed in Fig. 4b representing a decrease of 10% in the solid phase. It was used to represent the contour of the edematous tissue in the performed simulation. In addition, the area limited by this isoline and the domain boundary are used to quantify the edematous tissue and compare it to the patient edema area identified by the T2 mapping exam (Fig. 1b). The area of the edema measure by cardiac MRI is equal to 1.755cm2, whereas the area obtained by the numeric simulation is 1.728cm2. As one can observe, the relative error is equal to 1.56%.

Fig. 4
figure 4

Comparison between the patient edema, identified by a medical doctor, and the results of the numerical simulation. Panel a show the T2 mapping imaging exam, Fig. 1b, segmented into a binary image. Panel b show the fluid phase distribution after 10 h of simulation and a contour line around the simulated edema

Additional experiments were performed to verify if the random distribution of the lymph vessels affects the final result. Dozens of random distinct distributions of lymph vessels were generated using distinct seeds values, and it was observed that the impacts in the final result were insignificant.

Finally, once cbp is the parameter which couples the immune system to the hydro-mechanical model, it is important to quantify how the proposed model responds to different values of cbp. So, five simulations were performed varying the parameter cbp, ranging its value from 40.0 to 80.0, and the resulting edematous tissue area was evaluated for each simulation. The results of this analysis is reported in Table 5. As one can observe, a variation of ±20 in cbp values results in an edema 4.76% smaller for cbp=40, and 10,67% bigger for cbp=80.

Table 5 Influence of distinct cbp values on the edema area


Figures 3a and b present the numerical results of the simulation of interactions between the human immune system, represented by leukocytes, and pathogens. It is worth to notice that in the simulated scenario it is considered that pathogens start the infection on the epicardium (i.e. outside the domain boundary), more specifically in the boundary of the target edema showed in the T2 mapping exam (Fig. 1b). This infection spreads over time along the simulated myocardium. The simulation starts with no leukocytes in the tissue and after pathogens start to reproduce and spread, leukocytes arrive at the infected tissue from the bloodstream. Moreover, a scenario in which leokocytes and pathogens coexist was imposed to analyze the mechanical deformation of the simulated tissue, but it must be stressed that this immune model is able to reproduce different scenarios [51]. Furthermore, it is important to highlight that the parameter cbp is the core of the coupling strategy, once it measures the influence of the pathogen to the hydraulic permeability of the microvascular bed wall. So, since the fluid source in the modeled system is the capillary bed, the variation of edema size is directly affected by cbp as shown in Table 5.

Figures 3c and 3d present the results of the hydro-mechanical model, which simulate the mechanical deformation of the extracellular matrix and tissue cells, both defined as solid phase under the assumptions of a poroelastic media approach. It is important to notice that the fluid accumulation in the infection site is only possible due to the mechanical model, which in this case is triggered by the inflammatory response of the immune system model [49].

Using two types of cardiac MRI, T2 mapping and LGE, the size of the edema was quantitatively measured and highlighted in the MRI T2 mapping (Fig. 1b). Additionally, the imaging exam was segmented into a binary image, presented in Fig. 4a, representing the edematous area in grey and the remaining tissue in black, with the purpose of making easier the comparison between the clinical exam results and the computational ones. Figure 4b shows the fluid phase distribution after 10 h of simulation and a contour line around the simulated edema. The area of the simulated edema matches the area obtained by the T2 mapping exam with an error of 1.56%. One can observe that the position and shape of the simulated edema (Fig. 4b) are very similar to the clinical ones (Fig. 4a), which allow us to conclude that the proposed mathematical model is a very promising tool to model the edema formation in myocarditis.


This study presents a patient-specific model for edema formation in myocarditis. Using patient data obtained with cardiac MRI exams, we were able to create a mesh specific to fit his/her heart. The set of proposed PDEs, that couples the hydro-mechanical model with a mathematical description of the immune system reaction to an invading pathogen, was then solved using this mesh as the domain. The results were able to reproduce the size, position, and shape of myocardial edema. In particular, the edema area obtained by the solution of the numerical method was compared to the edema area obtained clinically.

The major contribution of this study is the quantitative validation of the new hydro-mechanical model using for this purpose clinical data. The obtained result is very promising, indicating that the mathematical model could reproduce some key aspects of edema formation. The quantitative comparison was possible due to the use of the T2 mapping cardiac MRI, a recent technique of cardiac exam which allows quantitative measures of the edematous area. We hope that, in the near future, models like the one presented in this study can be sufficiently precise for reducing the use of animal experiments in early stages of drug testing, as well as to improve drug safety. Additionally, we hope that a better understanding of the processes involved in edema formation in myocarditis can speed up the development of new treatments for patients that are diagnosed with this disease.

The present work could be extended in many distinct ways. The immune system model could be extended to include more details or to model a specific pathogen. In addition, the edema formation in other living tissues could be simulated. The mechanical model presented in this work has some limitations. For example, the influence of the fiber direction in the deformation is not taken into account. Moreover, it is also important the notice that the two-dimensional approach can not consider anisotropic advection-diffusion in the orthogonal direction of the plane axis. To overcome these limitations, in the near future we plan to extend the model to three dimensions using the whole patient-specific heart geometry. Furthermore, in this study we consider that the hydro-mechanical model is coupled with the infection model by means of Lp and σ functions, but we do not consider the feedback of the mechanical deformation influencing the infection model, i.e. we consider a one-way coupling strategy. In the near future, we also plan to analyze the influence of the mechanical feedback to the infection model.


Patient data and image processing

Cardiac MRI was performed in a 27-year-old male patient (85kg, 171cm 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.2mmol of gadolinium-based contrast using a sequence with these parameters: gradient-echo inversion recovery (GRE-IR) sequence (TR 6.0ms; TE 4.18ms; matrix 256×128; FOV 350×260; resolution 1.37×2.0×8mm; FA 25o; BW 130Hz/px; TI 400ms; 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. $$

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, Db is the diffusion coefficient of the pathogen through the interstitial fluid, qb denotes the pathogen reproduction, and rb 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 qb and rb 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} $$
$$\begin{array}{*{20}l} r_{b} &= \lambda_{nb} C_{l} C_{p}, \end{array} $$

where cp is the pathogen growth rate in the interstitial tissue; Cl 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}} $$

where \(C_{l}:\Omega \times I\rightarrow \mathbb {R}^{+}\) is the leukocytes concentration in the interstitial fluid, Dn is the diffusion coefficient of the leukocyte through the interstitial tissue, χnb is the leukocyte chemotaxis rate, qn is the leukocyte source, which represents the leukocyte extravasation from bloodstream to interstitium, and rn 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} $$
$$\begin{array}{*{20}l} r_{n} &= \lambda_{bn} C_{l} C_{p} + \mu_{n}C_{l}, \end{array} $$

where Cn, 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. $$

where qc and ql represent the capillary bed network and the lymphatic system influence in Ω, respectively, and q=(qc+ql)ρ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 ql 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), $$

where q0 is the normal lymph flow, P0 is the normal interstitial fluid pressure. Moreover, Vmax is the maximum lymph flow, Km 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 qc is given by:

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

where cf is the filtering coefficient given by Lp(S/V), Lp is the hydraulic permeability of the microvascular bed wall, and (S/V) is the vessel surface area per unit volume; P and Pc 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}), $$

where Lp0 is the health tissue hydraulic permeability of the capillary wall, Cp is the pathogen concentration, and cbp 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})}, $$

where σ0 is the oncotic reflection coefficient in a non-inflamed tissue, Cp is the pathogen concentration, and cbr 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. $$

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. $$

where vD 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. $$

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}. $$

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}), $$

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}, $$

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. $$

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. $$

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. $$

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).

Fig. 5
figure 5

a Finite element mesh generated based on Fig. 1b. This mesh represents a slice of the long axis of the left ventricle of the heart of the patient. b Distribution of triangular finite elements representing the lymph vessels (white) in the domain. The lymph vessels were uniformly and randomly placed over the domain corresponding to a total of 2.9% of the elements

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.

Availability of data and materials

The data sets supporting the results of this article are included within the article and the references. Model parameters and initial conditions used in simulations are included in this published article. Code to solve the mathematical model can be made available upon request to the authors.



Cardiac magnetic resonance


Finite element method


Human immune system


Interstitial fluid pressure


International society and federation of cardiology


Magnetic resonance imaging


Ordinary differential equation


Partial differential equation


Streamline upwind Petrov-Galerkin


United States


World Health Organization


  1. Richardson P, McKenna W, Bristow M, Maisch B, Mautner B, O’Connell J, Olsen E, Thiene G, Goodwin J, Gyarfas I, Martin I, Nordet P. Report of the 1995 world health organization/international society and federation of cardiology. task force on the definition and classification of cardiomyopathies. Circulation. 1996; 93:841–2.

    Article  CAS  PubMed  Google Scholar 

  2. Kindermann I, Barth C, Mahfoud F, Ukena C, Lenski M, Yilmaz A, Klingel K, Kandolf R, Sechtem U, Cooper LT, et al.Update on myocarditis. J Am Coll Cardiol. 2012; 59(9):779–92.

    Article  PubMed  Google Scholar 

  3. Basso C, Calabrese F, Corrado D, Thiene G. Postmortem diagnosis in sudden cardiac death victims: macroscopic, microscopic and molecular findings. Cardiovasc Res. 2001; 50(2):290–300.

    Article  CAS  PubMed  Google Scholar 

  4. Friedrich MG, Sechtem U, Schulz-Menger J, Holmvang G, Alakija P, Cooper LT, White JA, Abdel-Aty H, Gutberlet M, Prasad S, Aletras A, Laissy J-P, Paterson I, Filipchuk NG, Kumar A, Pauschinger M, Liu P. Cardiovascular magnetic resonance in myocarditis: A jacc white paper. J Am Coll Cardiol. 2009; 53(17):1475–87.

    Article  PubMed  PubMed Central  Google Scholar 

  5. Caforio AL, Pankuweit S, Arbustini E, Basso C, Gimeno-Blanes J, Felix SB, Fu M, Heliö T, Heymans S, Jahns R, et al.Current state of knowledge on aetiology, diagnosis, management, and therapy of myocarditis: a position statement of the european society of cardiology working group on myocardial and pericardial diseases. Eur Heart J. 2013; 34(33):2636–48.

    Article  PubMed  Google Scholar 

  6. Goldsby RA, Kindt TJ, Kuby J, Osborne BA. Immunology, 5th edn.New York: W. H. Freeman; 2002.

    Google Scholar 

  7. Abbas AK, Lichtman AH. Basic Immunology Updated Edition: Functions and Disorders of the Immune System, 4th edn.Philadelphia: Elsevier Health Sciences; 2012.

    Google Scholar 

  8. Sompayrac L. How the Immune System Works. Oxford: Wiley-Blackwell; 2012.

    Google Scholar 

  9. Staub NC. Pulmonary edema. Physiol Rev. 1974; 54(3):678–811.

    Article  CAS  PubMed  Google Scholar 

  10. Thompson RB, Pagano JJ, Chow K, Sekowski V, Paterson I, Ezekowitz J, Anderson T, Dyck JR, Haykowsky MJ, et al.Subclinical pulmonary edema is associated with reduced exercise capacity in hfpef and hfref. J Am Coll Cardiol. 2017; 70(14):1827–8.

    Article  PubMed  Google Scholar 

  11. Rockson SG. Lymphedema. Am J Med. 2001; 110(4):288–95.

    Article  CAS  PubMed  Google Scholar 

  12. Cornely ME. Lipedema and lymphatic edema. Berlin: Springer; 2006, pp. 10–14.

    Book  Google Scholar 

  13. Brunelle CL, Swaroop MN, Skolny MN, Asdourian MS, Sayegh HE, Taghian AG. Hand edema in patients at risk of breast cancer–related lymphedema: Health professionals should take notice. Phys Ther. 2018; 98(6):510–17.

    Article  PubMed  PubMed Central  Google Scholar 

  14. Laine G, Allen S. Left ventricular myocardial edema. lymph flow, interstitial fibrosis, and cardiac function. Circ Res. 1991; 68(6):1713–21.

    Article  CAS  PubMed  Google Scholar 

  15. Eitel I, Friedrich MG. T2-weighted cardiovascular magnetic resonance in acute cardiac disease. J Cardiovasc Magn Reson. 2011; 13(1):13.

    Article  PubMed  PubMed Central  Google Scholar 

  16. Puntmann VO, Zeiher AM, Nagel E. T1 and t2 mapping in myocarditis: seeing beyond the horizon of lake louise criteria and histopathology. Expert Rev Cardiovasc Ther. 2018; 16(5):319–30.

    Article  CAS  PubMed  Google Scholar 

  17. Kim PK, Hong YJ, Im DJ, Suh YJ, Park CH, Kim JY, Chang S, Lee H-J, Hur J, Kim YJ, et al.Myocardial t1 and t2 mapping: techniques and clinical applications. Korean J Radiol. 2017; 18(1):113–31.

    Article  PubMed  PubMed Central  Google Scholar 

  18. Spieker M, Katsianos E, Gastl M, Behm P, Horn P, Jacoby C, Schnackenburg B, Reinecke P, Kelm M, Westenfeld R, et al.T2 mapping cardiovascular magnetic resonance identifies the presence of myocardial inflammation in patients with dilated cardiomyopathy as compared to endomyocardial biopsy. Eur Heart J Cardiovasc Imaging. 2017; 19(5):574–82.

    Article  Google Scholar 

  19. Spieker M, Haberkorn S, Gastl M, Behm P, Katsianos S, Horn P, Jacoby C, Schnackenburg B, Reinecke P, Kelm M, et al.Abnormal t2 mapping cardiovascular magnetic resonance correlates with adverse clinical outcome in patients with suspected acute myocarditis. J Cardiovasc Magn Reson. 2017; 19(1):38.

    Article  PubMed  PubMed Central  Google Scholar 

  20. Siettos CI, Russo L. Mathematical modeling of infectious disease dynamics. Virulence. 2013; 4(4):295–306.

    Article  PubMed  PubMed Central  Google Scholar 

  21. Best K, Perelson AS. Mathematical modeling of within-host zika virus dynamics. Immunol Rev. 2018; 285(1):81–96.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Conover T, Hlavacek AM, Migliavacca F, Kung E, Dorfman A, Figliola RS, Hsia T-Y, Taylor A, Khambadkone S, Schievano S, et al.An interactive simulation tool for patient-specific clinical decision support in single-ventricle physiology. J Thorac Cardiovasc Surg. 2018; 155(2):712–21.

    Article  PubMed  Google Scholar 

  23. Oliveira RS, Alonso S, Campos FO, Rocha BM, Fernandes JF, Kuehne T, dos Santos RW. Ectopic beats arise from micro-reentries near infarct regions in simulations of a patient-specific heart model. Sci Rep. 2018; 8(1):16392.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  24. Bonin CRB, Fernandes GC, dos Santos RW, Lobosco M. A qualitatively validated mathematical-computational model of the immune response to the yellow fever vaccine. BMC Immunol. 2018; 19(1):15.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  25. Guerra A, Belinha J, Jorge RN. Modelling skin wound healing angiogenesis: a review. J Theor Biol. 2018.

    Article  CAS  PubMed  Google Scholar 

  26. Pearson JE, Krapivsky P, Perelson AS. Stochastic theory of early viral infection: continuous versus burst production of virions. PLoS Comput Biol. 2011; 7(2):1001058.

    Article  CAS  Google Scholar 

  27. Xavier MP, Bonin CR, dos Santos RW, Lobosco M. On the use of gillespie stochastic simulation algorithm in a model of the human immune system response to the yellow fever vaccine. In: 2017 IEEE International Conference on Bioinformatics and Biomedicine (BIBM): 2017. p. 1476–82.

  28. Beauchemin C, Samuel J, Tuszynski J. A simple cellular automaton model for influenza a viral infections. J Theor Biol. 2005; 232(2):223–34.

    Article  CAS  PubMed  Google Scholar 

  29. Bernaschi M, Castiglione F. Selection of escape mutants from immune recognition during hiv infection. Immunol Cell Biol. 2002; 80(3):307.

    Article  CAS  PubMed  Google Scholar 

  30. Kohler B, Puzone R, Seiden PE, Celada F. A systematic approach to vaccine complexity using an automaton model of the cellular and humoral immune system: I. viral characteristics and polarized responses. Vaccine. 2000; 19(7):862–76.

    Article  CAS  PubMed  Google Scholar 

  31. Bauer AL, Beauchemin CA, Perelson AS. Agent-based modeling of host–pathogen systems: the successes and challenges. Inf Sci. 2009; 179(10):1379–89.

    Article  Google Scholar 

  32. Celada F, Seiden PE. A computer model of cellular interactions in the immune system. Immunol Today. 1992; 13(2):56–62.

    Article  CAS  PubMed  Google Scholar 

  33. Chang ST, Linderman JJ, Kirschner DE. Multiple mechanisms allow mycobacterium tuberculosis to continuously inhibit mhc class ii-mediated antigen presentation by macrophages. Proc Natl Acad Sci U S A. 2005; 102(12):4530–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Vodovotz Y, Chow CC, Bartels J, Lagoa C, Prince JM, Levy RM, Kumar R, Day J, Rubin J, Constantine G, et al.In silico models of acute inflammation in animals. Shock. 2006; 26(3):235–44.

    Article  CAS  PubMed  Google Scholar 

  35. Su B, Zhou W, Dorman K, Jones D. Mathematical modelling of immune response in tissues. Comput Math Methods Med. 2009; 10(1):9–38.

    Article  Google Scholar 

  36. Flegg JA, Byrne HM, Flegg MB, McElwain DS. Wound healing angiogenesis: the clinical implications of a simple mathematical model. J Theor Biol. 2012; 300:309–16.

    Article  PubMed  Google Scholar 

  37. Pigozzo AB, Macedo GC, Dos Santos RW, Lobosco M. On the computational modeling of the innate immune system. BMC Bioinformatics. 2013; 14(6):7.

    Article  Google Scholar 

  38. Pigozzo AB, Missiakas D, Alonso S, dos Santos RW, Lobosco M. Development of a computational model of abscess formation. Front Microbiol. 2018; 9:1355.

  39. Quintela deBM, dos Santos RW, Lobosco M. On the coupling of two models of the human immune response to an antigen. BioMed Res Int. 2014; 2014.

    Article  CAS  Google Scholar 

  40. Cattaneo L, Zunino P. A computational model of drug delivery through microcirculation to compare different tumor treatments. Int J Numer Methods Biomed Eng. 2014; 30(11):1347–71.

    Article  CAS  Google Scholar 

  41. Jain RK, Martin JD, Stylianopoulos T. The role of mechanical forces in tumor growth and therapy. Ann Rev Biomed Eng. 2014; 16:321.

    Article  CAS  Google Scholar 

  42. Phipps C, Kohandel M. Mathematical model of the effect of interstitial fluid pressure on angiogenic behavior in solid tumors. Comput Math Methods Med. 2011; 2011.

    Article  Google Scholar 

  43. Reis RF, dos Santos RW, de Oliveira Campos J, Lobosco M. Interstitial pressure dynamics due to bacterial infection. Mecánica Computacional Bioeng Biomech (B). 2016; 34(1):1181–94.

    Google Scholar 

  44. Reis RF, dos Santos RW, Lobosco M. A plasma flow model in the interstitial tissue due to bacterial infection. Lect Notes Comput Sci. 2016:335–45.

    Chapter  Google Scholar 

  45. Biot MA. General theory of three-dimensional consolidation. J Appl Phys. 1941; 12(2):155–64.

    Article  Google Scholar 

  46. Selvadurai A, Suvorov A. Coupled hydro-mechanical effects in a poro-hyperelastic material. J Mech Phys Solids. 2016; 91:311–33.

    Article  Google Scholar 

  47. Suvorov A, Selvadurai A. On poro-hyperelastic shear. J Mech Phys Solids. 2016; 96:445–59.

    Article  Google Scholar 

  48. Berger L, Bordas R, Burrowes K, Grau V, Tavener S, Kay D. A poroelastic model coupled to a fluid network with applications in lung modelling. Int J Numer Methods Biomed Eng. 2016; 32(1).

    Google Scholar 

  49. Reis RF, Weber dos Santos R, Rocha BM, Lobosco M. On the mathematical modeling of inflammatory edema formation. Comput Math Appl (accepted for publication in March 2019)E. 2019; 0(0):0–0.

    Google Scholar 

  50. Reis RF, Rocha BM, dos Santos RW, Lobosco M. IEEE International Conference on Bioinformatics and Biomedicine, BIBM 2018, Madrid, Spain, December 3-6, 2018 In: Zheng HJ, Callejas Z, Griol D, Wang H, Hu X, Schmidt H, Baumbach J, Dickerson J, Zhang L, editors.. IEEE Computer Society: 2018. p. 1418–24.

  51. Pigozzo AB, Macedo GC, Weber dos Santos R, Lobosco M. Computational modeling of microabscess formation. Comput Math Methods Med. 2012; 2012.

    Article  Google Scholar 

  52. Kellman P, Aletras AH, Mancini C, McVeigh ER, Arai AE. T2-prepared ssfp improves diagnostic confidence in edema imaging in acute myocardial infarction compared to turbo spin echo. Magn Reson Med. 2007; 57(5):891–7.

    Article  PubMed  PubMed Central  Google Scholar 

  53. Salemi VMC, Rochitte CE, Shiozaki AA, Andrade JM, Parga JR, de Ávila LF, Benvenuti LA, Cestari IN, Picard MH, Kim RJ, Mady C. Late gadolinium enhancement magnetic resonance imaging in the diagnosis and prognosis of endomyocardial fibrosis patients. Circ Cardiovasc Imaging. 2011; 4(3):304–11.

    Article  PubMed  Google Scholar 

  54. Guyton AC, Hall JE. Textbook of Medical Physiology, 12 edn. Guyton Physiology Series. Philadelphia: Elsevier Saunders; 2006.

    Google Scholar 

  55. Keener JP, Sneyd J. Mathematical Physiology vol. 8. New York: Springer; 1998.

    Book  Google Scholar 

  56. Starling EH. On the absorption of fluids from the connective tissue spaces. J Physiol. 1896; 19(4):312–26.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Scallan J, Huxley VH, Korthuis RJ, Vol. 2. Capillary Fluid Exchange: Regulation, Functions, and Pathology. United States of America: Morgan & Claypool Life Sciences; 2010, pp. 1–94.

    Google Scholar 

  58. Cheng A. H-D. Poroelasticity vol. 27. Switzerland: Springer; 2016.

    Book  Google Scholar 

  59. Logg A, Mardal K-A, Wells GN. Automated Solution of Differential Equations by the Finite Element Method. Lysaker: Springer; 2012.

    Book  Google Scholar 

  60. Langtangen HP, Logg A. Solving PDEs in Python: The FEniCS Tutorial I vol. 1. Oslo: Springer; 2017.

    Google Scholar 

  61. Donea J, Huerta A. Finite Element Methods for Flow Problems, 1st edn.West Sussex: John Wiley & Sons; 2003.

    Book  Google Scholar 

  62. Geuzaine C, Remacle J-F. Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities. Int J Numer Methods Eng. 2009; 79(11):1309–31.

    Article  Google Scholar 

  63. Reis RF, dos Santos RW, Lobosco M. Influence of the immune system on the biological dynamics of the interstitial fluid pressure. Singapore: Springer; 2017. pp. 304–7.

    Chapter  Google Scholar 

  64. Rahier J-F, De Beauce S, Dubuquoy L, Erdual E, Colombel J-F, Jouret-Mourin A, Geboes K, Desreumaux P. Increased lymphatic vessel density and lymphangiogenesis in inflammatory bowel disease. Aliment Pharmacol Ther. 2011; 34(5):533–43.

    Article  PubMed  Google Scholar 

Download references


Not applicable.

About this supplement

This article has been published as part of BMC Bioinformatics Volume 20 Supplement 6, 2019: Towards computational modeling on immune system function. The full contents of the supplement are available online at


The authors would like to thank the infrastructure provided by Universidade Federal de Juiz de Fora (UFJF) and Instituto de Ensino e Pesquisa José Michel Kalaf (IEPJMK) to perform this study. This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) - Brasil - Finance Code 001 (scholarship), Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) - Brasil (scholarship) and Fundação de Amparo à Pesquisa do Estado de Minas Gerais (FAPEMIG) - Brasil (equipments). This manuscript has not received sponsorship for publication costs.

Author information

Authors and Affiliations



Conception and design of the mathematical model: RFR, BMR, RWS and ML. Computational implementation of the mathematical model: RFR. RFR has also been involved in drafting the manuscript. JLF was responsible for the analysis and interpretation of exams. JLF and TRS participated in writing, providing perspective from the Cardiology standpoint. All authors read, revised and approved the final manuscript.

Corresponding author

Correspondence to Ruy Freitas Reis.

Ethics declarations

Ethics approval and consent to participate

The collection of images and written informed consent for this study was undertaken under a larger project that was approved by the local Ethics Committee (Comitê de Ética em Pesquisa Hospital Municipal Dr. Mario Gatti - CEP 711/2007, Certificado de Apresentação para Apreciação Ética - CAAE 0509.0.146.000-07) for investigating patients with different etiologies for cardiomyopathies.

Consent for publication

All patient data were anonymized before evaluation. All patients filled an institutional form consenting to participate of a study investigating patients with different etiologies for cardiomyopathies.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver( applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Reis, R., Fernandes, J., Schmal, T. et al. A personalized computational model of edema formation in myocarditis based on long-axis biventricular MRI images. BMC Bioinformatics 20 (Suppl 6), 532 (2019).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI:


  • Computational immunology
  • Myocarditis
  • Poroelasticity
  • Mathematical modeling
  • Biomechanics