An Introductory Overview of Image-Based Computational Modeling in Personalized Cardiovascular Medicine

Cardiovascular diseases account for the number one cause of deaths in the world. Part of the reason for such grim statistics is our limited understanding of the underlying mechanisms causing these devastating pathologies, which is made difficult by the invasiveness of the procedures associated with their diagnosis (e.g., inserting catheters into the coronal artery to measure blood flow to the heart). Likewise, it is also difficult to design and test assistive devices without implanting them in vivo. However, with the recent advancements made in biomedical scanning technologies and computer simulations, image-based modeling (IBM) has arisen as the next logical step in the evolution of non-invasive patient-specific cardiovascular medicine. Yet, due to its novelty, it is still relatively unknown outside of the niche field. Therefore, the goal of this manuscript is to review the current state-of-the-art and the limitations of the methods used in this area of research, as well as their applications to personalized cardiovascular investigations and treatments. Specifically, the modeling of three different physics – electrophysiology, biomechanics and hemodynamics – used in the cardiovascular IBM is discussed in the context of the physiology that each one of them describes and the mechanisms of the underlying cardiac diseases that they can provide insight into. Only the “bare-bones” of the modeling approaches are discussed in order to make this introductory material more accessible to an outside observer. Additionally, the imaging methods, the aspects of the unique cardiac anatomy derived from them, and their relation to the modeling algorithms are reviewed. Finally, conclusions are drawn about the future evolution of these methods and their potential toward revolutionizing the non-invasive diagnosis, virtual design of treatments/assistive devices, and increasing our understanding of these lethal cardiovascular diseases.


INTRODUCTION: IMAGE-BASED MODELING OF THE HEART
Heart disease is the leading cause of death in the U.S., with one person dying from it every 37 s, or about 647,000 each year (i.e., 1 in every 4 deaths), and amounting to a $219 billion per year burden to the public health system (CDC, 2019). Understanding it is very difficult, because it is a complex interaction of biomechanics, electrophysiology and non-Newtonian hemodynamics. This is further complicated by the interaction with external medical devices (pacemakers, pumps, etc.) that are commonly implanted in order to assist a failing or dysfunctional heart. Moreover, the heart's properties (e.g., shape, structure, stiffness, electrical conductivity) that play an important role in determining its pumping ability are patient specific. Finally, it is difficult to extract information about the physiological processes occurring in living hearts, due to its constant motion, and the fact that invasive probing can be life threatening. For these reasons, Image-Based Modeling (IBM) -a patient-specific experimentally constrained computational approach -is a lucrative way for gaining novel insight into the cardiovascular diseases and their treatments.
An illustrative example of the IBM's usefulness is HeartFlow Inc. -a company located in California, United States with about 300 employees and backed by $467 million capital investment (Craft, 2019). Their application is the diagnosis of Coronary artery disease (CAD) -an impairment of blood flow in the arteries that supply the heart, due to cholesterol plaque buildup. The disease is one of the most misdiagnosed: a recent study, which included data from more than 1,100 U.S. hospitals, found that over half of the more than 385,000 patients with suspected CAD underwent an invasive coronary angiography (ICA) only to find out that they did not have the disease (Patel et al., 2014). This is bad because ICA is an invasive technique that in itself could lead to mortality, because it uses catheters inserted into the femoral (groin) or radial (wrist) arteries to measure pressure difference across a coronary artery stenosis in order to check the likelihood of a blockage's presence.
Conversely, HeartFlow calculates pressure differences virtually by simulating the blood flow through the patients' own arteries, the structure of which is derived from a threedimensional (3D) computerized tomography (CT) scan. This information is then used to calculate the fractional flow reserve (FFR), which is a statistic used to assess the hemodynamic significance of the stenosis by determining the ratio of the pressures before and after the narrowing. Therefore, this technology effectively serves as a non-invasive alternative to the ICA (Figure 1) (HeartFlow, 2019).
The HeartFlow's method has been evaluated in four large prospective clinical trials, enrolling a total of more than 1,100 patients at major medical centers worldwide. It received the European Economic Area CE mark in 2011 and U.S. FDA clearance in August 2019 (i.e., it is currently commercially available in the U.S.) (FDA, 2019). To date, clinicians have used the HeartFlow approach for over 30,000 patients in the diagnosis of heart disease (Kim, 2019). Therefore, it serves as the most mature IBM application in the context of cardiovascular disease. Yet, it is also one of the simplest in that it does not include the heart itself, and the blood assumed a homogeneous (i.e., no cells) fluid. At the same time, more advanced models are coming online as well. Yet, they are relatively unknown outside of this niche field.
To that end, this review is organized as follows: Section "Methods: Literature Search" describes our literature search methods; Section "Background: Cardio Electromechanics and Hemodynamics" provides a brief background of the relevant cardiovascular physiology that explains how the tissue electromechanics and blood biology are interrelated in vivo; Section "Geometry Module" reviews how the model geometry is obtained using imaging, in order to establish a connection with the individual's unique anatomy; Sections "Electrophysiology Module, " "Biomechanics Module, " "Simplified Hemodynamics Module, " and "Hemodynamics with Thrombogenesis Module" illustrate the mathematical formulation of the simulation modules used by these models, such as electrophysiology, biomechanics, simplified hemodynamics with and without thrombogenesis, respectively. Finally, Section "Summary and Conclusion" presents our summary and conclusions regarding the inputs and outputs of all cardiac modules as well as the directions that the field of personalized cardiovascular IBM is expected to go into.

METHODS: LITERATURE SEARCH
In order to provide a "big picture" snapshot overview of the image-based cardiovascular modeling and its potential penetration into the clinical sector, we gathered works from the recent (i.e., approximately the last 5 years) proceedings of the following meetings: Interagency Modeling and Analysis Group Consortium Meetings, National Institute for Mathematical and Biological Synthesis workshops, Personalized Medicine Coalition resources, International Workshop on Cancer Systems Biology meetings and conferences, Pacific Symposium on Biocomputing, Biomedical Engineering Society, American Institute of Chemical Engineers. Additionally, we performed manual searches with key words and terms including "image-based modeling of the heart, " "patient-specific cardiovascular modeling, " "cardio electrophysiology/biomechanics/electromechanics/hemodynamics image-based modeling" or "ventricular thrombosis modeling, " etc. for both research and review articles on databases, such as PubMed central, Web of Science, Research Gate, and Google Scholar. Additionally, we relied on the corresponding author's own decade and a half long experience of working on biomedical image-based simulations. The obtained publication database was then screened by running citation reports to identify groups of researchers (typically led by a senior professor, who is joined by collaborators, postdocs, and students) that have established a track-record of being active within the various sub-areas of the cardiovascular modeling fields. Furthermore, to avoid bias (and to keep the work manageable) we tried to limit the literature sampling to just one most relevant publication from each of the groups. However, this was not always possible, because some of the researchers dominate their respective niches; and have published more than one article critical to our review.

Macroscopic Overview of How the Three Physics Are Coupled With Each Other
Before going into the details of the computational models, it is first important to understand the three types of coupled physics occurring in the heart: electrical signal conduction, biomechanics of the contraction and hemodynamics (which could also include clot formation and embolism). Figure 2A illustrates the cardiac conduction system (CCS)a heterogeneous complex 3D network of highly specialized conductive cells (SA node, AV node, bundle of His, bundle branches, and Purkinje fibers) that transfer signals through the heart and cause it to contract. The electrical activity is initiated at the sinoatrial node (i.e., the natural pacemaker of the heart) where voltage signals called "action potentials" are produced periodically. Next the signals travel to the AV node, through the Bundle of HIS, down its branches and through the Purkinje Fibers. Ultimately, they are propagated to the myocardium (i.e., the muscular tissue of the heart) through discrete sites called Purkinje-Myocyte Junctions (not shown), causing the left and the right ventricles to contract independently of each other. This creates a double pumping action of the blood (see Figure 2B). Specifically, oxygen-poor blood returns from the body to the right side of the heart (i.e., atrium and ventricle), which then sends it to the lungs for re-oxygenation. Oxygen-rich blood from the lungs then enters the left side of the heart and is pumped through the aorta back to the body. The blood can also carry thrombi (or their embolized pieces) from other parts of the body, and/or the clots could be generated within the heart itself via activation of platelets (blood cells responsible for clot formation) and the coagulation cascade (a series of biochemical reactions that results in the formation of a polymer mesh that facilitates the structural integrity of the clot). The presence of these objects in the cardiovascular system can interfere with the mechanical action of the heart by creating rigid obstructions. Furthermore, the blood clots can also get stuck in the cardio-vasculature and block the delivery of metabolites to the heart tissue. This leads to necrosis of the latter, commonly referred to as an ischemia or a heart attack.

Structural Importance of the Myocardium
The organization of the cardiomyocyte fibers in the heart's walls is thought to be critical to both the conductive and to the mechanical properties of the organ. Specifically, the contractile myocytes cells that cause the heart ventricles to beat are arranged in fibers (see inset of Figure 2B). These fibers make up the walls of the heart, which are lined with collagen and elastin extracellular matrix on the inside (i.e., endocardium) and the outside (i.e., epicardium). Their thickness varies both spatially and temporarily throughout the cardiac cycle.
If one were to take a representative sample from the left ventricle (which pumps the most blood) as in Figure 3A, it would be possible to see that the 3D layered organization of the myocytes changes throughout the wall thickness from the epicardium to the endocardium. In fact, Figure 3A shows that the muscle fiber direction rotates from +50 • to +70 • (sub-epicardial region) to nearly 0 • in the mid-wall region to −50 • to −70 • (sub-endocardial region) with respect to the circumferential direction of the left ventricle (Holzapfel and Ogden, 2009). Finally, Figure 3B show that the myocyte fibers (or myofibrils) are arranged into composite layers (or sheets), which are interconnected by collagen fibers. Therefore, for IBM to be physiologically representative, it must account for how this intricate structure affects the complex physics that occur in the heart.

GEOMETRY MODULE
The most common ways for obtaining a realistic macroscopic morphology of the heart and its surrounding blood vessels are computerized tomography (CT) and magnetic resonance imaging (MRI). However, the typical MRI/CT machines provide relatively coarse resolution datasets of the personalized cardiac geometry, with large gaps between slices (Schulte et al., 2001;Frangi et al., 2002;Appleton et al., 2005). This necessitates the use of interpolation procedures. Hence, the microscopic details (e.g., blood vessels, trabeculations, the Purkinje Fibers, the location and activity of the PMJs, the orientations of the myofibril sheets, etc.) are harder to resolve due to their small size. Yet, they strongly determine the electrophysiological and biomechanical properties of cardiac tissue (Watson et al., 2018). Consequently, there are three main methods for accounting for these fine details:

Rule-Based Heuristics
The most rudimentary approach is to generate these features mathematically, based on observed trends (see Figure 4A) (Streeter Daniel et al., 1969). Briefly, the longitudinal fiber direction is assumed to rotate clockwise from the endocardium to the epicardium. Specifically, it is made parallel to the long axis of the papillary muscles, trabeculae at these regions and parallel to the endocardial and epicardial surfaces at the ventricular  walls. Lastly, the fiber orientation in the septum is assumed to be running along the ventricular walls (Bayer et al., 2012). A popular way to personalize the algorithm to a patient specific structure of the heart is to use the minimal distance between the imaged endocardial and the epicardial surfaces to approximate orientation of the fibers. More stable and advanced rule based approaches, such as the Laplace-Dirichlet method also exist (Bayer et al., 2012). However, the heuristics are not guaranteed to be physiologically accurate, nor are they fully patient specific. Yet they remain the most common approach due to their low cost and ease of implementation, as well as due to the difficulty of imaging the microstructural details in a beating heart in vivo.
And, as Figure 4 shows, they yield results that are comparable with the best of the imaging techniques (which typically require for the heart to be explanted and fixed in order to acquire such fine details).

Histology/Optical Microscopy
The fiber orientation can also be approximated from histology of explanted hearts (Vetter and McCulloch, 1998;Deng et al., 2012), where the tissue is sliced into very thin 2D sections and dyed using special agents that highlight the features of interest. Although it is possible to create a 3D reconstruction based on the 2D slices using this approach, manual sectioning of the tissue may result in uneven slice thicknesses and feature distortion (Burton et al., 2006). Additionally, confocal microscopy has been used to image the fine structure of the myocardium (Hooks et al., 2002(Hooks et al., , 2006. However, the light penetration depth into the sample is typically limited to ∼100 mm. Hence, imaging a whole heart using this technique is also impractical. Therefore, both histology and confocal are typically used to provide localized information on explanted samples only. However, this information is useful for validating the rule-based methods, the in vivo imaging, and the modeling results.

Diffusion-Tensor MRI, Micro-CT With Contrast and 3D Ultrasound Backscatter Tensor Imaging
Additional detail can be obtained from MRI images using a Gadolinium contrast agent (Bishop et al., 2010) and a special technique called Diffusion Tensor imaging (DTI) (see Figure 4B). The latter maps the diffusion of water molecules in the biological tissues, which is not free, but reflects the interactions with obstacles like macromolecules, fibers and membranes. For the cardiac DTI, it is well known that the direction of the primary eigenvector corresponding to each voxel of the received images matches the longitudinal axis of cardiac myocytes (Scollan et al., 1998;Holmes et al., 2000). This information can then be mapped onto the volumetric mesh of a 3D cardiac macrogeometry to include the microscopic fiber orientation (Plank et al., 2009;Vadakkumpadan et al., 2010). Likewise, microcomputed tomography (mCT) with iodine staining has also recently been used to assess the myocyte fiber orientation in the heart tissue (Aslanidi et al., 2013). However, both techniques are too slow to capture a beating heart in 3D without motion artifacts. Luckily, advanced ultrasound-based imaging techniques are coming online, which can map the myocardial fibers orientation and its dynamics with a temporal resolution of 10 ms during a single cardiac cycle, non-invasively and in-vivo in entire volumes (Papadacci et al., 2017). However, given the novelty, complexity and cost of these techniques, they are not yet widely available to the majority of the cardiovascular IBM researchers.

Imaging-to-Modeling Pipeline
Perhaps the most difficult aspect of the in vivo scanning of live hearts is the need to perform significant image alignment using "registration" techniques. Furthermore, once the images are aligned, they must be "segmented" to identify the various tissue types and structural features of interest within the data. The segmentation can be either based on contrast dyes and/or on morphological feature detection (both manual and automated). The images can also be enhanced using digital post-processing, such as deconvolution (i.e., minimizing noise caused by objects outside of the imaging plane) and structure tensor analysis (e.g., enhancing visibility of the structural features for the fiber orientation detection) (Burton et al., 2006;Zhao et al., 2013). Numerical interpolation and machine learning techniques can further enhance the apparent resolution of the images digitally. Finally, the cardiovascular geometry must be "meshed" (i.e., broken up into pieces) to discretize the objects obtained from the images as a set of finite elements for numerical analysis. The latter is a simulation necessity that enables solving systems of complex (e.g., partial differential) equations by recasting them as algebraic approximations of the true solution. The trade-off for the simplified math is that the solution accuracy must be increased by making the mesh finer. This grows the number of equations that must be solved simultaneously, and thus the computational resource and time requirements. Overall, the imaging pipeline procedures are often very complicated and necessitate manual labor. This is both cumbersome (e.g., due to the large size of the high-resolution images) and subjective (e.g., due to the lack of contrast agents which necessitate user input). Therefore, there is an on-going effort to automate the imaging-to-modeling pipeline (Bishop et al., 2010).

ELECTROPHYSIOLOGY MODULE
The simplest types of the heart models tend to be focused on cardiac arrhythmias. This is an "umbrella" term for irregularities in the conduction or pacing of the electrical signals that control the heartbeat rhythm. Given that some arrhythmias can lead to mortality, it is important to understand the underlying electrophysiological mechanisms of these disorders. Yet, electrocardiograms of the heart provide only limited information, which often fails to predict lethal outcomes (Goldberger et al., 2011). Therefore, computational modeling offers a better alternative for studying these diseases.
Most arrhythmia models focus on the electrophysiology of the heart, while assuming that it is isolated from the biomechanics of the contractions that lead to the pumping of the blood. As mentioned in Section "Methods: Literature Search, " the contraction of cardiomyocytes is initiated by electrical impulses called "action potentials, " which travel through the cardiac conduction system (see Figure 2A) into the myocardium. This electrical potential travels from one cell to another in the form of ions that pass through gap junctions between the cardiomyocytes (see Figure 5).
The cardiomyocytes are polarized, meaning that there is an electrical potential across the cell membrane: in the resting state the cells are more negative on the inside and positive on the outside, while the charge polarity is temporarily reversed as the action potential passes through them. This reversal occurs via the transport of Ca ++ and Na + ions from the outside of the cells to their inside, and K + ions in the reverse direction (see Figure 5). The internalization of the Ca ++ ion is especially important to the contraction of the cardiomyocytes, because it triggers a sub-cellular signaling cascade that generates tension inside of the cells. Therefore, the electrophysiology models simulate the propagation of the ionic currents through the myocardium.
However, given that there are many cells in the heart, and each one of them has the ionic channels and transmembrane potentials, the problem is an inherently multiscale one (see Figure 6). On the subcellular scale (see Figure 6-LEFT), differential equations are used to model the transport of ionic species based on the Hodgkin-Huxley formulation (originally developed for the propagation of action potentials in neurons) (Hodgkin and Huxley, 1952). The subcellular models are then combined into cell scale models that can account for up to dozens of different ionic species and signaling intra-cellular cascades (see Figure 6, CENTER). Among these, the leading model is currently considered to be by O'Hara et al. (2011), which is based on experimental data from >150 undiseased human hearts. Finally, the individual ion currents are used to calculate the overall transmembrane potential, and its transport across the myocardium is treated as a diffusion across a homogenous medium (i.e., no discrete cells are considered by these models).
There are two types of formulations for the organ-level diffusion (which is related to the ionic conductivity) of the action potential across the myocardium: (1) the bidomain formulation, which considers different diffusivities inside and outside of the cell (Quarteroni et al., 2017) (see Figure 6, RIGHT): The bidomain formulation is used for simulating the action potential propagation throughout the myocardium in the intra-and the extra-domains separately. The myocardium is assumed to be a continuum in which the potential is considered to vary along the longitudinal direction of the conducting cells, while it is constant in the transversal (or radial) directions (Quarteroni et al., 2017). In this formulation, v i , v e and 'v' are intracellular, extracellular and transmembrane potentials, respectively; i i app and i e app stand for applied stimuli on the intraand extracellular spaces, respectively; i ion are the ionic currents following a Hodgkin-Huxley-type description for different ionic species (Hodgkin and Huxley, 1952); w are gating variables taking values in [0,1] that regulate the transmembrane currents and have a mutual relation with the intracellular concentrations c of different ionic species (which also vary depending on the values of transmembrane potentials v) (Quarteroni et al., 2017); C m is the membrane capacitance; 'χ' is the ratio of membrane area per tissue volume; D i and D e are the conductivity tensors of the intraand extracellular media, respectively. and (2) the monodomain formulation, which simplifies the problem by considering only the transmembrane potential (Quarteroni et al., 2017): In this formulation, the cardiac tissue is also assumed to be a continuum, but the current conservation is written in terms of the transmembrane potential v only (i.e., not considering the intraand extracellular potentials) (Quarteroni et al., 2017). Instead, the intracellular and extracellular diffusivities are assumed to be proportional to each other, and therefore can be represented by a single variable. Herein, D 0 is the conductivity tensor in a fixed reference state, J is the determinant of the deformation gradient tensor, which represents the volume change of a deformable object. The trade-off for the simplicity is that the monodomain model is unable to describe cardiomyocyte repolarization patterns. For this reason, the bidomain model is more widely used (Bishop and Plank, 2011). Table 1 summarizes the most common electrophysiology module personalization approaches encountered in the recent IBM works, while Figure 7 maps the relationships between the module's inputs, outputs and applications. In this, and in the subsequent modules, the heart's macroscopic anatomy can be personalized for a specific patient by acquiring its geometry from the in vivo imaging (see "Geometry Module" section). Furthermore, pathological tissue remodeling (e.g., locations and extension of infarct scars, diffuse fibrosis, etc.) can be accounted for via the imaging as well (Vadakkumpadan et al., 2009;Mewton et al., 2011;Dass et al., 2012;Ashikaga et al., 2013;Arevalo et al., 2016;Trayanova et al., 2017). However, there are two important types of microscopic structural information that cannot be completely personalized yet: the myocardial fiber orientation and the CCS. Both are, unfortunately, still too difficult to image in a moving heart in vivo. Thus, typically an approximated personalization is performed to include these microscopic features using the rule-based algorithms (see "Geometry Module" section). Additional CCS approximation methods are based on: early activation points obtained from the literature (Durrer et al., 1970); manual delineation of CCS on the endocardial surfaces (Romero et al., 2010); ex vivo data obtained by means of histological studies of animal hearts ; from in vivo electro-anatomical maps (EAMs) (Cardenes et al., 2014;Palamara et al., 2014). The EAMs data can provide the location of some of the PMJs, and can also be used to reconstruct patient-specific electrical activation patterns (Cardenes et al., 2014;Palamara et al., 2014).

Module Personalization
Likewise, the conductivity values of different tissue zones (normal and abnormal/damaged) are typically too difficult to measure in living human patients, and are thus initialized based on accepted literature values: 2-3 m/s in the His-Purkinje system and 0.3-0.4 m/s in the conductive myocardial cells (Ideker et al., 2009). They are then either further adjusted to match the human myocardium conduction velocity (CV) (i.e., speed at which action potentials are distributed throughout the tissue) measured using explanted hearts Trayanova et al., 2017;Lopez-Perez et al., 2019), or are tuned to fit patient-specific electrical activation patterns obtained from electrocardiograms (ECG), body surface potential maps (BSPM) or EAM (Sermesant et al., 2008;Lopez-Perez et al., 2015). Unfortunately, at the cellular level, the patient-specific transmembrane current dynamics (i.e., i ion ) cannot yet be measured; and hence, the existing mathematical models are not personalized at such detail. Similarly, the electrical heterogeneity between the different regions (e.g., transmural heterogeneity in the ventricular walls), the electrical remodeling and the effects due to an individual's genetic mutations on the cardiac electrophysiology cannot yet be accounted for (Lopez-Perez et al., 2015). However, a cellular level electrophysiology model that best matches the patient's pathology can instead be chosen from either existing literature datasets that are representative of a patientgroup (Krueger et al., 2013a,b) or from patch-clamp studies Geometry with orientation and position of heart Geometry was personalized and tissue conductivity values were adjusted to match patient-specific ECG data. Kayvanpour et al., 2015 MRI (with torso geometries) Geometry with orientation and position of heart Geometry was personalized and tissue conductivity values were adjusted to match the CV in human cardiac tissue. Geometry was personalized and tissue conductivity values were adjusted to match patient-specific activation mapping data using EAMs.

Roney et al., 2018
Cardiac magnetic resonance procedure Detailed geometry, including the atrial structure, aortic arch, caval veins, torso surface, trabeculated myocardium between the wall and the intracavitary blood.
Geometry was personalized and tissue conductivity values were tuned to match patient-specific activation mapping using EAMs and ECG data. Potse et al., 2014 of cells harvested from pathologic zones of the patient (Cabo and Boyden, 2003;Decker and Rudy, 2010). Additionally, the extracellular ion concentrations can be estimated and set into a model from personalized measurements of blood electrolyte concentrations (Krueger et al., 2013a,b).

Module Outputs and Applications
Overall, the electrophysiology modeling studies the normal conduction in the heart, as well as the pathological mechanisms that arise and cause cardiac arrhythmias. It is typically used to calculate physiological parameters (see Figure 7) like: the Mean firing rate (i.e., the number of spikes during a cardiac cycle divided by cycle duration, spike/s) (Behradfar et al., 2014); Re-entrant arrhythmias propagation Deng et al., 2016;Trayanova and Chang, 2016;Prakosa et al., 2018) (i.e., a propagation of an impulse that fails to die out after normal activation of the heart and continues to re-excite it after the refractory period has ended, Antzelevitch, 2001); Phase singularities Pathmanathan and Gray, 2018;Roney et al., 2018) which represent the sites in which the activation state cannot be determined, because the particular location is surrounded by activation states ranging from fully activated to fully recovered (Valderrabano et al., 2003); Activation rate gradient which quantifies how fast the transmembrane voltage V m changes in different cardiac regions (Smith et al., 2011;Boyle et al., 2013;Potse et al., 2014); Successful retrograde propagation which measures whether conduction at a terminal Purkinje node is successful or refractory (Trayanova, 2011;Behradfar et al., 2014); and the organization of electrical wavelets as they propagate through the myocardium (Starobin et al., 1996;Keldermann et al., 2009;Trayanova, 2014). In our experience, the majority of the electrophysiological modeling is used for elucidating the mechanisms of cardiac arrhythmia, especially for "reentrant propagation of complex waves" (e.g., effects of cardiac microstructure, spiral wave breakup, early afterdepolarizations, scroll-wave filaments, action potential duration, electrical alternans, etc.) (Trayanova, 2011); as well as for prediction of arrhythmia risks in specific patients. Furthermore, these models are also used to examine the mechanisms of defibrillation shock in the heart for terminating arrhythmia, as well as for increasing the understanding of ablation targets in the arrhythmia treatments (Trayanova et al., 2017).

Module Personalization Example
The following is a discussion of a representative example of the personalized electrophysiology modeling applied to an arrythmia risk assessment in post-infarcted hearts. Specifically, personalized 3D computer models of the post-infarction hearts was constructed based on clinical MRI of specific patients. First, an individualized geometric model of the postinfarction ventricles was reconstructed from late-gadolinium-enhanced-MRI , with representations of both the scar and the infarcted border zones. Due to the difficulty imaging the myocardial fiber orientation from a moving heart in vivo, an approximated personalization was performed using a rule-based algorithm (Bayer et al., 2012). Region-specific cell and tissue electrical properties were then assigned to the electrophysiological model based on literature data. After that, a virtual multi-site delivery of electrical stimuli from various biventricular locations was conducted, in order to computationally determine all the ventricular tachycardia reentrant pathways that the infarct-remodeled ventricular substrate can sustain. This methodology was then validated in an arrhythmia risk prediction clinical study including 41 patients and significantly surpassed several existing clinical metrics in predicting upcoming arrhythmic events ).

BIOMECHANICS MODULE
The next level of complexity are the models of myocardial abnormalities/heart failure (HF) and the blood pumping assist devices such as the Left-Ventricular Assist Device (LVAD). Since these models are interested in how the presence of abnormalities or assist devices affects the blood circulation, they must account for the contraction solid mechanics and the blood hemodynamics. However, if they are not interested in clot formation, the models are simplified by homogenizing the blood flow. Therefore, they are not as complicated as the ones that do account for the thrombosis. Yet, they are still complex, because they include solving multiple different types of coupled physics (see Figure 8, LEFT) applied in different parts of the heart (see Figure 8, RIGHT). Since the electrophysiology in these models is treated similarly to the arrhythmia models, the cardio biomechanics framework will be discussed next. Central to this physics type is the fact that the myocytes in the heart contain rod-like structures called myofibrils, which are composed of repeating contractile units called "sarcomeres" (see Figure 9, LEFT). Each sarcomere contains thin and thick filaments, made from actin and myosin proteins, respectively.
After depolarization, calcium enters the cardio-myocytes through the ion channels in their membrane, and triggers the release of cytosolic calcium stored in the sarcoplasmic reticulum (a storage compartment) of the cells through a cascade of intracellular signaling (see Figure 9, RIGHT). This release of the internal calcium leads to the binding of myosin heads to the actin filaments in the sarcomeres, which in turn causes the filaments to slide against each other and contract the entire cell (see bottom inset in Figure 9, LEFT). This process is called the "crossbridge mechanism" and it corresponds to the active tension in the biomechanical calculations of the heart.
The kinetics of this process have been modeled using Monte Carlo (Walcott and Sun, 2009) and partial differential equations (Huxley, 1957). However, both approaches are too computationally expensive to be calculated at the organ level. Consequently, the latter model has been simplified by considering a single cross bridge representative of the whole distribution (i.e., mean field theory) (Negroni and Lascano, 2008;Rice et al., 2008;Washio et al., 2012) or by averaging the distributions over a single cell (Bestel et al., 2001). Simpler yet is the assumption that the active contraction of the individual cells depends on the intracellular ionic concentrations and the local deformation gradient (Quarteroni et al., 2017). Ultimately, the microscopic sarcomere sliding velocity is related to the macroscopic strain along the myocardial fibers via a constitutive relationship (i.e., microscopic rate-of-strain depends on the macroscopic strain), such as the Hill-Maxwell rheological model which is a modification from Hill's force-velocity relation (Sainte-Marie et al., 2006). Herein, a specific choice of the attachment and detachment rates was modified so that they were not only dependent upon the sarcomere strain, but also on the strain rate.
The active contraction can be expressed by the following active stress formulation from Rossi et al. (2014): And the active strain formulation can be derived from Quarteroni et al. (2017): Where, W A is the active component of the free energy, F A is the active deformation, 'c' is the intracellular calcium concentration, I 4,f is the local deformation gradient invariant in the myofiber direction, ∂I E 1 and ∂I E 4,f are elastic invariants (described in more details in Rossi et al., 2014). The incorporation of the microscopic  active tension,F A , at the tissue level is a key aspect in the multiscale framework of the cardiac function (Quarteroni et al., 2017) and can be derived by: where R F−L (I 4,f ) is a function that represents the force-length relationship of the cardiac cells (defined in Ruiz-Baier et al., 2014) which depends on I 4,f ; f (c) specifies the amount of force generated by the cross-bridges in response to intracellular calcium release; and α is a positive parameter. Thus, theF A represents the active tension generated within the sarcomeres, which then drives the macroscopic muscular contractions. In addition to the active tension generated by the cross-bridge mechanism, the solid mechanics calculations must also account for the passive stiffness of the myocardium (e.g., when it is expanded by the blood flow entering the heart). This is modeled by assuming that the myocardium is an isotropic linear elastic tissue. The general formulation of passive stress (i.e., Cauchy stress) can be expressed by Avazmohammadi et al. (2019): Where, the second Piola-Kirchhoff stress tensor 'S' can be described in terms of the stored energy density function 'W' through: S = 2∂W/∂C, and 'W' can be derived based on Holzapfel-Ogden constitutive law which is divided into three parts -the isotropic isochoric part, the isotropic volumetric part, and the orthotropic part: where J = det(F), C = F T F, I 1 = tr(C),Ī 1 = J −2/3 I 1 , I 4,f = Cf 0 · f 0 ,Ī 4,f = J −2/3 I 4,f , I 4,s = Cs 0 · s 0 , I 4,s = J −2/3 I 4,s , and the material parameters a, a f , a s , a fs , b, b f , b s , b fs are experimentally fitted (Quarteroni et al., 2017). The parameter 'κ' is the bulk modulus that "penalizes" local volume changes to enforce the incompressibility of the tissue. Furthermore, since stretching a cell changes the distance between the gap junctions and their neighbors, this leads to changes in the ion channels, and consequently, in the conductivity of the action potential from cell to cell. This electromechanics coupling is typically included by modifying the conductivity tensor in the original equation of the electrophysiological propagation (see Equation 2). Specifically, the fixed reference state conductivity tensor 'D 0 ' is replaced with the spatial configuration conductivity tensor 'D'. Furthermore, an explicit dependence on the solid deformation tensor 'F' is included into the conductivity tensor in order to account for the geometric feedback, due to deformation of the tissue structure (Quarteroni et al., 2017): Moreover, an additional inward ionic current, induced by the stretch-activated channels, also contributes to the depolarization. This is commonly modeled as follows: Where, 'E' and 'g' are the reversal potential and the maximal conductance of the channels. Lastly, in vivo the heart is immersed in a pericardial fluid; and it is also loosely supported by a flexible double-layered pericardium membrane. Therefore, spring-like external support enforced by Robin-type boundary conditions (Moireau et al., 2012) are typically tuned to mimic the global motion of the modeled heart. An explicit solution of the pericardium-heart contact problem has been proposed as well (Fritz et al., 2014). Table 2 summarizes the most common biomechanics module personalization approaches encountered in the recent IBM publications, while Figure 10 the relationships between the module's inputs, outputs and applications. Before personalizing the modeling of the passive (i.e., resting) properties of the myocardium, the parameters of the Holzapfel and Ogden model (see Equations 7 and 8) are typically initialized using experimental data from biaxial (Krishnamurthy et al., 2013) and shear tests (Dokos et al., 2002;Sommer et al., 2015) of explanted myocardial tissue. The passive mechanics parameters are then further optimized to match the patient-specific end-diastolic pressure and volume (EDPV) relations (Krishnamurthy et al., 2013;Meoli et al., 2015;Finsberg et al., 2018;Palit et al., 2018). The EDPV relation is a graphical representation of the pressurevolume loop related to the passive filling of the left ventricle during diastole (i.e., relaxation), and is a measure of the passive ventricle stiffness. The chamber volume of the left ventricle at the end of the diastole is defined as the end-diastolic volume (EDV), which is used to estimate the preloading volume of the heart and indicate the stiffness of the ventricle.

Module Personalization
The active contraction modeling, on the other hand, is coupled with the cardiac hemodynamics, so it is important to link it with the patient-specific hemodynamic metrics corresponding with the end-systolic condition. For example, the pumping ability of the heart is represented by the end-systolic pressure (ESP) and the end-systolic volume (ESV), which are the peak values of pressure and ventricular volume at the end of systole (i.e., contraction), respectively. The active contraction is also associated with ejection fraction (EF), which is defined as a measurement of the percentage of blood leaving out of the left ventricle with each contraction. Therefore, the optimization procedure is performed to find the patient-specific parameters that best replicate the clinical hemodynamics of the patient in terms of the mean, maximum and minimum values of the pressures, flows and the cardiac volumes (such as the clinically recorded values of the ESP, ESV, and EF).

Module Outputs and Applications
Overall, the biomechanics models tend to perform stress analysis, which is used to evaluate the effects that a defective myocardium structure has on the heart function, and design new therapies and treatment devices for reducing the abnormal stress (Guccione et al., 2003;Wall et al., 2006;Lee et al., 2013;Finsberg et al., 2018) (see Figure 10). Additionally, they calculate the stiffness of the heart, which strongly corresponds to its ability to function normally and can be used as an indicator for HF (e.g., heart attacks caused by a diastolic dysfunction that occurs when the ventricle becomes too stiff or weak to pump blood effectively). Ultimately, however, since the goal of these models is to obtain the relationship between the biomechanics of the myocardium tissue and the blood pumping ability of the heart, they also include a simplified (i.e., no discrete blood cells) hemodynamics module, which is discussed in the next section.

Module Personalization Example
The following is a discussion of a representative example of the personalized Biomechanics module applied to five HF failure patients from San Diego Veteran's Affairs Medical Center. Specifically, personalized 3D models of ventricular biomechanics in their failing hearts were derived from cardiac CT imaging. The human fiber orientation was modeled using DT-MRI data from an isolated (i.e., fixed) human organ-donor heart, and then transposed to the specific patient's geometric model. The biomechanics model was then developed for optimizing the passive material properties to match previous experimental results on cardiac tissues and patient-specific end-diastolic pressure and volume relations. The material properties of the active contraction were also optimized to match patientspecific measured peak left ventricular pressures and end-systolic volumes. These components were then integrated to generate a multi-scale computational approach for the patient-specific hearts. The simulation results in the patients demonstrated good agreement with their measured echocardiographic and cardiac output parameters, such as EF and peak cavity pressures. This model was developed for stress analysis in HF patients and could be further developed with the goal of predicting treatments for heart disease under different interventions.

SIMPLIFIED HEMODYNAMICS MODULE
In the cardiovascular models where clot formation is not considered, the blood flow is simulated using the incompressible Navier-Stokes equations. This means they do not account for discrete cells floating in the plasma. Instead the blood is treated as a homogeneous weakly non-Newtonian fluid (Shibeshi and Collins, 2005), which flows mostly in the laminar regime [though the strong vortices can create transition to turbulence with Reynolds numbers in the 1500-2500 range (Quarteroni et al., 2017)]. The fluid-structure interaction between the blood and the myocardium walls is typically modeled explicitly using moving mesh approaches: for example, Arbitrary Lagrangian-Eulerian (Cheng et al., 2005;Chnafa et al., 2014;Su et al., 2014), immersed boundary (Kohl et al., 2001;Vigmond et al., 2008) and levelset based methods (Mihalef et al., 2011). The heart valves, on the other hand, are commonly approximated using the Bernoulli equation for orifice flow (Flachskampf et al., 1990;Vandervoort Pieter et al., 1995;Donati et al., 2017).
Additionally, the myocardium hemodynamics are typically coupled to the rest of the body's circulation via the Windkessel circuit model, which mimics the arterial blood pressure's waveform. This is a relatively simple method used to obtain the relationship between the blood flow and the pressure in a modeled segment through the resistive 'R' and the capacitance 'C' properties of the arterial vasculature (see Figure 11) (Wall et al., 2006). Specifically, the heart and the systemic arterial system are treated as a closed hydraulic circuit, which, contains a pump connected to a chamber partially filled with a liquid. As it is pumped, the latter compresses the air pocket in the chamber, which in turn pushes the liquid back out (i.e., creating a back-and-forth cycle). Consequently, the arterial compliance, the peripheral resistance, and the inertia are modeled as a capacitor, a resistor, and an inductor, respectively. In this model, the physiological variables such as pressure 'P' and flow 'Q' only vary as a function of time 't' (Morris et al., 2016). The relationship between the flow rate Q(t) and the pressure P(t) can be expressed by: It essentially states that the volumetric flowrate must equal to the sum of the volume stored in the capacitive element and the volumetric outflow through the resistive element. During the diastole there is no blood inflow (Q = 0), so the Windkessel equation can be solved for the P(t): where t d is the time of the start of the diastole and P(t d ) is the blood pressure at that time. Due to its simplicity, the Windkessel equation is frequently used to approximate various components and boundary conditions in the cardiovascular system (Morris et al., 2016). However, this model is only a rough approximation Geometry was personalized; passive mechanics parameters were fitted to match the previous experimental results on cardiac tissues and match patient-specific EDPV relations; active mechanics parameters were adjusted to match the patient-specific measured peak left ventricular pressures and ESV. Krishnamurthy et al., 2013 DT-MRI from explanted heart Fiber direction Cine cardiac MRI Geometry with ventricular dimensions and ejection fraction (EF).
Geometry was personalized; passive mechanics parameters were fitted to match the previous experimental results and match the previous results of EDPV relations (due to the unavailability of subject-specific ventricular pressures that require invasive measurements).

Palit et al., 2018
Cine cardiac MRI Geometry with ventricular dimensions and regional strain-time Geometry was personalized; passive mechanics parameters were fitted to match patient-specific EDPV relations; active mechanics parameters were adjusted to match the patient-specific end-systolic state.

Finsberg et al., 2018
MRI flow tracing and echocardiographic Doppler velocity tracings Geometry with ventricular dimensions and blood flow velocities.
Geometry was personalized; passive mechanics parameters were fitted to match the previous experimental results on cardiac tissues and match patient-specific EDPV relations; active mechanics parameters were adjusted to match the patient-specific measured peak left ventricular pressures. Meoli et al., 2015 MRI Geometry with ventricular dimensions Geometry was personalized; passive parameters were tuned so that the resultant EDV matched the corresponding MRI-measured cavity volume; active parameters were adjusted so that the resultant ESV matched the corresponding MRI-measured cavity volume. Lee et al., 2013 Cine MRI and Echocardiography Geometry with ventricular dimensions and EF Geometry was personalized; passive and active parameters were determined by minimizing the sum of the squared differences between computed and measured EF, stroke volume, EDV, ESV, end-diastolic pressure and end-systolic pressure. Kayvanpour et al., 2015 Cine MRI Geometry with ventricular volume waveforms Geometry was personalized; biomechanics parameters were tuned along with hemodynamics parameters; passive parameters were adjusted to match the measured EDPV; active parameters were adjusted to match patient-specific volume and pressure waveforms.

HEMODYNAMICS WITH THROMBOGENESIS MODULE
The last class of the cardiovascular models is the one in which the hydrodynamics of the discrete blood cells, and the biochemistry associated with their physiological activity, are of interest.
For example, such simulations could be focused on studying pathological clot formation inside of the heart (Choi et al., 2015; FIGURE 11 | Two-element Windkessel circuit analogy illustrated. Mittal et al., 2016;Seo et al., 2016;Harfi et al., 2017), or embolism into it from other parts of the body. Therefore, these types of models must simulate the blood as a suspension of deformable cells (e.g., platelets and/or red blood cells), whose fate is intertwined with the mechanical motion of the myocardium (and by association with the electrophysiology of the heart). This is typically done using Stokesian Dynamics methods, Dissipative Particle Dynamics, Completed Double Layer Boundary Integration Equation Method and Lattice Boltzmann Method (Wang and King, 2012). These are mesoscopic offand on-lattice frameworks that calculate trajectories of the cells under the influence of hydrodynamic and Brownian forces; while the deformation of the structure is typically simulated using continuum-based models that treat the cell membrane and intracellular fluids as homogeneous materials: some popular approaches are the Boundary Integral Method, the Immersed Boundary Method, and the Fictitious Domain Method (Li et al., 2017).
To make matters even more complicated, the mechanism of the blood clot formation strongly depends on the following three processes: Receptor-Ligand Binding, Platelet Activation and the Coagulation Cascade. Specifically, the initiation of thrombus development starts with tethering of circulating platelets onto the exposed subendothelial layer where a blood vessel is injured. This process involves bonding between the various receptors on the platelet surfaces to the extravascular proteins, such as the von Willebrand Factor. It is typically modeled using a Monte Carlo approach called Adhesive Dynamics, while the rate of the receptor-ligand bond formation and breakage are determined by the Bell Model that calculates the probability of dissociation events occurring over a specific timespan.
Once the platelets have been recruited to the injury site, they undergo a metamorphosis that is generally described as "activation." During this process, the membrane receptors transmit signals to the inside of the cells, which results in the dumping of chemical agonists stored in the internal vesicles called lysosomes and granules (e.g., dense and alpha). The release of these molecules then activates other neighboring platelets, which ultimately leads them to becoming "stickier" and compacting into the blood clot's body. Unfortunately, the bottom-up description of the process contains numerous unknowns (Dolan and Diamond, 2014). For this reason, it has been instead described via a top-down Neural Network approach, which was trained on patient-specific experimental data (Flamm et al., 2012).
Lastly, the coagulation cascade is a system of coupled biochemical reactions with two different initiation pathways, which both ultimately lead to the polymerization of soluble fibrinogen (a blood plasma protein) into an insoluble fibrin mesh that holds the clot together. The kinetic portion of the cascade involves 34 differential equations, with 42 rate constants, that cumulatively account for 27 independent equilibrium expressions and fates of 34 chemical species (Hockin et al., 2002). Additionally, the mass transport of these species must be tracked under the flow conditions experienced in the cardiovascular system (which involves Knudsen diffusion within the porous clot). Figure 12 summarizes the coupling of the various hemodynamics submodules, as well as the methods used to solve them. Table 3 summarizes the most common hemodynamics module personalization approaches encountered in the recent IBM works, while Figure 13 maps the relationships between the module's inputs, outputs and applications. Similarly to the previous modules, the macroscopic geometry of the heart, and that of the surrounding blood vessels, is commonly obtained for use in the Hemodynamics modules via the in vivo imaging techniques (such as MRI and CT). Additionally, for the body's circulation, the parameters values for the Windkessel model (i.e. 'R' and 'C' elements) are typically chosen to match the patient specific cardiac output, flow waveforms and pressure pulses (Kim et al., 2009(Kim et al., , 2010Kung et al., 2014) obtained via contrast-enhanced CT scans, Doppler ultrasound scanning and invasive blood pressure measurements (IBPM). Particularly, the Doppler ultrasonography allows the measurement of the cardiac output and heart rate; and the IBPM allows to measure the systolic and diastolic pressures (Bonfanti et al., 2019). Specifically, they are first chosen based on literature data (Segers et al., 2002;Kim et al., 2009) and are then iterated until the calculated FIGURE 12 | Overview of the specific methods, modules and models commonly used in the state-of-the-art multiscale platelet hemodynamics and thrombus development studies. (Adopted with permission from Wang and King, 2012). flow parameters match the subject's unique physiological profile (Kim et al., 2009).

Module Personalization
Unfortunately, due to the enormous complexities and computational costs of the heart modeling, oversimplifications are common in the personalization of the blood properties within the Hemodynamics modules. Given that it is typically assumed to be a Newtonian (or weakly non-Newtonian) fluid, it effectively does not contain discrete blood cells, such as the platelets or the erythrocytes. This, in turn, means that only simple flow properties, like fluid viscosity can be personalized to the subject. Consequently, the patient-specific biology of these cells (e.g., thrombotic propensity, sickle cell anemia, etc.) are omitted. Yet, in the non-cardio blood modeling, attempts at personalizing such functionality and disorders have been made: for example, the application of the Neural Networks trained on the patientspecific experimental data in order to phenotype the platelet activation (Flamm et al., 2012). However, we could not find an example of such an extensive blood personalization in the cardiovascular modeling literature. Figure 13 summarizes the inputs, outputs and applications of the cardiac hemodynamics module. Overall, the models that use the simplified hemodynamics formulation tend to calculate the blood flow parameters, like the pressure-volume relationship of the cardiac cycle. These, in turn, help to elucidate quantities that are key to understanding the heart dysfunctions: such as the end-diastolic and the end-systolic pressure volume relationships. Additionally, Fractional Flow Reserve, which is defined as the pressure difference across a coronary artery stenosis (e.g., a narrowing due to atherosclerosis), can be calculated to determine the likelihood that the latter would impede oxygen delivery to the heart and lead to a myocardial ischemia. Furthermore, the biomechanical models measure "compliance, " which is the ability of the blood vessel walls to stretch in order to accommodate an increasing amount of blood; and "resistance" (defined as the ratio of the pressure drop and the flow change across the segment) that the blood flow experiences due to viscous stresses and constrictions by the blood vessels. Most importantly, the biomechanics-hemodynamics modeling can be used to predict the left ventricular ejection fraction (LVEF) -a main indicator of HF, which is expressed as a percentage of how much blood the left ventricle pumps out with each contraction. Lastly, such models can be used to investigate the effects of the blood pumping assist devices (e.g., LVAD) on the cardiac function, which may otherwise be too difficult or costly to study experimentally. MRI Coronary arteries geometry Geometry was personalized; parameter values of lumped heart model (for the inlet) and Windkessel models (for the outlet) were adjusted to match subject-specific flow distribution and measured brachial artery pulse pressure. Kim et al., 2009 CT angiogram Right coronary artery with aneurysmal region Geometry was personalized; Windkessel model' parameters were adjusted to match subject-specific flow distribution and pressure at outlet boundary of the coronary artery. Kung et al., 2014 Phase contrast MRI 3-component flow velocity at two slice locations in the coronary aneurysm geometry

Module Outputs and Applications
Contrast-enhanced CT Aortic dissection geometry Geometry was personalized; inlet flowrate was obtained by adjusting a typical ascending aorta blood flow waveform to match the patient-specific hemodynamic data including cardiac output, heat-rate and systolic-to-diastolic duration ratio; for the outlet, Windkessel models parameters were adjusted to achieve patient-specific physiological flow distribution at each outlet, and to obtain the measured systolic and diastolic pressures at the inlet.

Bonfanti et al., 2019
Doppler ultrasonography Cardiac output and heart rate CT angiogram Coronary arteries geometry Geometry was personalized; parameter values of lumped heart model (for the inlet) and Windkessel models (for the outlet) were adjusted to match subject-specific flow distribution and measured aortic pressure.
Grande Gutierrez et al., 2017 4D cardiac CT images and echocardiogram Whole-heart geometry with kinematics of left ventricle lumen: heart rate, EDV, ESV, EF, stroke volume.
Geometry was personalized; the immersed boundary-based method is used to match the patient-specific kinematics of the left ventricular lumen including heart rate, EDV, ESV, EF and stroke volume. Harfi et al., 2017 Frontiers in Bioengineering and Biotechnology | www.frontiersin.org Conversely, the models that use the thrombogenesis hemodynamics module tend to try to assess the propensity of clot formation at (or near) the heart, based on parameters that are related to platelet activation: such as the Wall Shear Rate; the blood Residence Time Near a Damaged Tissue; the Ejection Fraction (i.e., the percentage of the blood leaving the left ventricle each time it contracts); Washout Ratio (i.e., the ratio of delayed ejection volume to the total ventricular blood at the beginning of the cycle). Additionally, some of these models simulate how drugs and clot breakup devices (e.g., Vena cava filters) help to protect the heart from pathogenic events. Overall, this is the most computationally expensive model type, due to the complexity of the thrombogenesis/embolism processes (which are themselves still being actively investigated) (Wang and King, 2012).

Module Personalization Example
The following is a discussion of a representative example of the personalized hemodynamics modeling applied to predicting the thrombosis risk in patients with Kawasaki disease (KD) (Grande Gutierrez et al., 2019). Thrombosis is a major complication associated with coronary artery aneurysms (CAAs) resulting from the KD. In this research, the aneurysm hemodynamics were investigated for thrombotic risk stratification in ten KD patients, and were compared to the standard clinical guidelines for anticoagulation therapy. The patient-specific models were generated from MRI data by performing an angiography of: the heart, the main coronary arteries (right, left main, left anterior descending and circumflex), and the aorta and its arch branches (the brachiocephalic artery, the left common carotid artery and the left subclavian artery). This was done via the injection of gadolinium-based contrast with a cardiac and respiratory-gated 3D TrueFISP sequence. Computational hemodynamic simulations were then performed in the reconstructed anatomical model using SimVascular software (Grande Gutierrez et al., 2019). The pulsatile flow, deformable arterial walls and Windkessel parameters were tuned to match the patient-specific arterial pressure and cardiac output. Local hemodynamics variables were derived from the simulation results, including the time-averaged wall shear stress, low wall shear stress exposure and blood residence time. These variables were then used to develop a framework for predicting the thrombosis risk. Although platelet activation and aggregation are typically associated with regions of higher fluid shear (Casa et al., 2015) and longer blood stagnation (Hathcock James, 2006), this study showed that a combination of low shear stress coupled with a high residence time correlated to thrombus formation in the KD CAAs patients. Furthermore, it was shown that the prediction of the thrombotic risk using the hemodynamic variables was validated with a higher sensitivity and specificity in comparison with the standard clinical metrics. In conclusion, this type of personalized computational modeling can be used to provide a non-invasive thrombotic risk stratification that is FIGURE 14 | Summary of cardiovascular IBM's modules and their respective inputs and outputs with corresponding references: 1 Krishnamurthy et al., 2013;Lopez-Perez et al., 2015;2 Krishnamurthy et al., 2013;Cardenes et al., 2014;Kayvanpour et al., 2015;Lopez-Perez et al., 2015;3 Krishnamurthy et al., 2013;Grande Gutierrez et al., 2017;Finsberg et al., 2018;4 Meoli et al., 2015;Bonfanti et al., 2019;5 Arevalo et al., 2016;Trayanova et al., 2017;Lopez-Perez et al., 2019;6 Krishnamurthy et al., 2013;Finsberg et al., 2018;Palit et al., 2018;7 Kung et al., 2014;Grande Gutierrez et al., 2017;Bonfanti et al., 2019;8 Choi et al., 2015;Seo et al., 2016;Harfi et al., 2017;9 Arevalo et al., 2016;Deng et al., 2016;Trayanova and Chang, 2016;Prakosa et al., 2018;10 Arevalo et al., 2016;Trayanova et al., 2017;Deng et al., 2019;11 Voorhees and Han, 2015;Walmsley et al., 2017;Lee et al., 2018;12 Koo et al., 2011;Min et al., 2012Min et al., , 2015Zhang et al., 2014;Min et al., 2015. more accurate than the current clinical approaches. This, in turn, can assist the long-term medical management of the KD patients with the CAAs.

SUMMARY AND CONCLUSION
Most of the published cardiovascular modeling reviews are typically oriented at an expert audience, which makes it difficult for the outsiders to understand the full medical potential of these methods. One of the barriers to penetrating the field is that these works tend to focus on just one or two specific aspects of the simulation approach at a time: such as imaging (Weese et al., 2013;Lamata et al., 2014;Watson et al., 2018), electrophysiology (Lopez-Perez et al., 2015;Rodriguez et al., 2015;Beheshti et al., 2016;Gray and Pathmanathan, 2018;Ni et al., 2018), biomechanics (Wang et al., 2015;Chabiniok et al., 2016), hemodynamics , electro-biomechanical coupling (Trayanova, 2011(Trayanova, , 2012Tobon-Gomez et al., 2013;Niederer et al., 2019b), biomechanicshemodynamics coupling (Tang et al., 2010;Sun et al., 2014), etc. In contrast, our manuscript provides a "big picture" overview of the components that these models are built from; their mechanisms, inputs, outputs and connecting pipelines; the underlying physiological processes that they represent; their image-based personalization to the individual patient's unique anatomy; and their applications to the different cardiovascular disease understanding and treatments.
As a part of our review, it was found that although this type of modeling holds a tremendous potential for revolutionizing personalized cardiovascular medicine, it is still in its infancy (with HeartFlow being the only commercially available product). Furthermore, due to the slow speed of high-resolution imaging, most of the IBM in academia still rely on scans of dead hearts (as opposed to beating ones). This, however, is expected to change, as the imaging speeds of mCT and MRI are increased. As far as the physics being modeled, the modules, their inputs and outputs are summarized in Figure 14. The simplest application of the cardiovascular IBM is to the study of arrythmias, which simulates the propagation of electrical impulses through the myocardium, while ignoring the biomechanics and hemodynamics. At the subcellular level, these models calculate transfer of ions across the cell membrane channels, while at the macro level the transfer of potential between the cells is modeled as a diffusive process. Conduction irregularities, ablation targets, and effects of defibrillation are just some of the outputs provided by the electrophysiological models.
The more complex models account for the biomechanical processes occurring in the myocardium. These models are focused on how abnormalities in the tissue morphology and stiffness affect the pumping efficiency of the heart. At the subcellular level, they simulate how the electrically driven calcium signaling governs the crossbridge mechanism of the active tension generation via the actin-myosin interactions in the cells' sarcomeres. Additionally, the solid mechanics calculation also involves the passive stress of the myocardium, which represents the stress-strain relationship of the cardiac fibers without the electrical stimulation. The resulting tissue deformation is then backwards-coupled to the electrophysiology, since the stretching of the cells can change the gap junction distance to their neighbors (which results in changes to the ionic conductivity).
Lastly, these models are coupled to the hemodynamics calculations via fluid-structure interactions. If only the blood flow without the clot formation is of interest, then the structure of the fluid is approximated to be homogeneous and Newtonian; while the rest of the body's circulation is introduced as an oscillating pressure boundary condition that is assumed to behave like a simple RC circuit. This type of model can provide insight into HF due to deformities and obstructions, as well as allow for virtual design and testing of assisted pumping devices. On the other hand, if clotting information is necessary, the blood must be treated as a suspension of deformable particles, with receptor ligand interactions, intra-cellular signaling, and coupled biochemical reactions representing the coagulation cascade. Although a lot more involved, this type of model can be useful for antithrombotic drug development and design of clot breakup devices (e.g., Vena cava filters) meant to protect the heart. However, for simplicity most cardiovascular IBM do not account for the hydrodynamics and biomechanics of the individual blood cells. Likewise, the coagulation cascade and the platelet activation, both of which are central to thrombogenesis, are oversimplified relative to the state-of-the-art within the non-cardiovascular blood modeling field.
Overall, the cardiovascular IBM is expected to become more mainstream as the computational and imaging technologies advance. This could potentially revolutionize how cardiovascular medicine is done in the future. Yet, significant improvements are still required to personalize the models more. For example, a common limitation across all of the modules is that the myocardium (e.g., conductivity and tissue stiffness) and the blood (e.g., hematocrit, coagulation cascade and platelet activation kinetics and deficiencies) properties that they use are typically estimated based on empirical measurements performed ex vivo and using samples that are not derived from the same individual (or even from the same species for that matter). Therefore, better imaging methods need to be developed, such that these properties could be estimated by scanning the individual patient. Furthermore, a finer image resolution is needed to capture the individualized variations in the conductive and contractile fibers, and their junctions. Likewise, the computational methods need to improve their modeling of the intra-cellular processes (e.g., the crossbridge mechanism) for the cardiovascular IBM to become more physiologically representative and adopted by the mainstream clinical market. However, given the fast pace of the technological progress, the near future impact of the IBM on the cardiovascular medicine is imminent.

AUTHOR CONTRIBUTIONS
TN and OK devised the manuscript, performed review of literature, and wrote the manuscript. RV provided the feedback and edits on all aspects of the manuscript. All the authors contributed to the article and approved the submitted version.

FUNDING
The study was financially supported by the New Jersey Health Foundation Grant #PC101-17 and the Gustavus and Louise Pfeiffer Research Foundation's Major Investment Grant.

Monodomain Formulation
A mathematical framework that is used in Electrophysiology modeling to simulate the propagation of the Action Potential through the myocardium. As opposed to the Bidomain Formulation (see definition), it does not consider both the intra-and the extra-cellular domains of the ion exchange across the cardiomyocyte membrane. Instead, the framework is expressed in terms of a transmembrane potential, which assumes that the intracellular and extracellular ion diffusivities are proportional to each other.

Neural Network
A set of artificial intelligence algorithms, modeled loosely after how the human brain learns, that are designed to recognize patterns within complex datasets. In cardiovascular modeling specifically, they are used to represent processes that are too difficult to describe mathematically using first-principle methods: for example, platelet "activation" -a cascade of cell membrane surface receptor activation, intracellular signaling, dumping of granules, etc. -in the context of thrombosis.

Stored (a.k.a., Strain) Energy Density Function
Energy per unit volume of a material temporarily deformed by an applied force, like a coiled spring or a stretched elastic band. This quantity is commonly used to formulate constitutive force-deformation relationships that characterize the spatial and temporal variations in the orthotropic (i.e., different in the axial, radial, and circumferential directions) properties of the myocardium. It is largely dominated by the structural arrangements of the myocardial fibers in the heart tissue.