In Silico Clinical Trials for Cardiovascular Disease

The rapid development of information technologies, simulation software packages, and medical devices in recent years provides the opportunity for collecting a large amount of clinical information. Creating comprehensive and detailed computational tools has, therefore, become essential to process specific information from the abundance of available data.

From the physicians' point of view, it is of paramount importance to distinguish "normal" versus "abnormal" phenotypes in a specific patient to estimate disease progression, therapeutic responses, and future risks. Recent computational models have significantly improved the integrative understanding of the behavior of heart muscles in hypertrophic (HCM) and dilated (DCM) cardiomyopathies¹. It is crucial to use a high-resolution, detailed, and anatomically accurate model of whole-heart electrical activity, which necessitates massive computation times, dedicated software, and supercomputers^1,2,3. A methodology for a real 3D heart model has been recently developed using a linear elastic and orthotropic material model based on Holzapfel experiments, which can accurately predict the electrical signal transport and displacement field within heart⁴. The development of novel integrative modeling approaches could be an effective tool for distinguishing the type and severity of symptoms in patients with multigenic disorders and assessing the degree of impairment in normal physical activity.

There are, however, many new challenges for patient-specific modeling. The physical and biological properties of the human heart are not possible to fully determine. Non-invasive measurements usually include noisy data from which it is difficult to estimate specific parameters for the individual patient. Large-scale computation requires a lot of time to run, whereas the clinical time frame is limited. Patient personal data should be managed in such a way that generated metadata can be reused without compromising patient confidentiality. Despite these challenges, multi-scale heart models can include a sufficient level of detail to achieve predictions that closely follow observed transient responses, thereby providing promise for prospective clinical applications.

However, regardless of the substantial scientific effort by multiple research labs and the significant amount of grant support, currently, there is only one commercially available software package for multiscale and whole heart simulations, called SIMULIA Living Heart Model⁵. It includes dynamic electro-mechanical simulation, refined heart geometry, a blood flow model, and complete cardiac tissue characterization, including passive and active characteristics, fibrous nature, and electrical pathways. This model is targeted for use in personalized medicine, but the active material characterization is based on a phenomenological model introduced by Guccione et al.^6,7. Therefore, SIMULIA cannot directly and accurately translate the changes in contractile protein functional characteristics observed in numerous cardiac diseases. These changes are caused by mutations and other abnormalities at the molecular and subcellular levels⁶. The limited use of SIMULIA software for a small number of applications in clinical practice is a great example of today's struggles in developing higher-level multiscale human heart models. On the other hand, it motivates the development of a new generation of multiscale program packages that can trace the effects of mutations from the molecular to organ scale.

The main aim of electrophysiology of the heart is to determine signal propagation inside the torso and the properties of all compartments^4,5,6. The SILICOFCM⁸ project predicts cardiomyopathy disease development using patient-specific biological, genetic, and clinical imaging data. It is achieved with multiscale modeling of the realistic sarcomeric system, the patient's genetic profile, muscle fiber direction, fluid-structure interaction, and electrophysiology coupling. The effects of left ventricle deformation, mitral valve motion, and complex hemodynamics give detailed functional behavior of the heart conditions in a specific patient.

This article demonstrates the use of the SILICOFCM platform for a parametric model of the left ventricle (LV) generated automatically from patient-specific ultrasound images using a fluid-structure heart model with electromechanical coupling. Apical view and M-mode view analyses of LV were generated with a deep learning algorithm. Then, using the mesh generator, the finite element model was built automatically to simulate different boundary conditions of the full cycle for LV contraction⁹. On this platform, users can directly visualize the simulation results such as pressure-volume, pressure-strain, and myocardial work-time diagrams, as well as animations of different fields such as displacements, pressures, velocity, and shear stresses. Input parameters from specific patients are geometry from ultrasound images, velocity profile in the input and output boundary flow conditions for LV, and specific drug therapy (e.g., entresto, digoxin, mavacamten, etc.).

The protocol in this study was approved by the UK National Health Service Health Research Authority North East-Tyne & Wear South Research Ethics Committee with the reference number 18/NE/0318 on 6 February 2019 and was adopted by the Institutional Review Board of each participating center. The study was conducted within the principles of Good Clinical Practice and following the Declaration of Helsinki. Informed consent was obtained from all subjects involved in the study. The patient information is kept anonymous.

1. Workflow for ultrasound M-mode or Apical view DICOM image analysis and LV parameter extraction

NOTE: To begin this protocol, the user must log in to the SILICOFCM platform⁸ and choose the appropriate workflow (i.e., ultrasound image analysis using M-mode or Apical view). The workflow for ultrasound M-mode DICOM image analysis and LV parameter extraction, for example, involves several steps. The first step is template matching. Templates to be matched include all the relevant necessary borders. It should be emphasized that, for every new dataset, the template should be extracted manually just once for each dataset that corresponds to a specific ultrasound machine. The workflow diagram is presented in Figure 1. The area matching the template will be extracted from the analyzed image after Felsenszwalb's efficient graph-based image segmentation algorithm finds "strong" boundaries that correspond to the borders of the septum and the wall of the LV. Based on these borders and the places where the heart diameters are the largest (corresponds to diameter in diastole) and the smallest (corresponds to diameter in systole), various LV dimensions will be calculated. The user should define whether the view is M-mode or Apical view.

1. Log in to the platform with a username and password. Under Virtual Population Module, choose ultrasound M-mode or Apical view workflow.
2. In the list of available workflows, select the ultrasound-merged workflow.
3. In the File Upload section, upload images and DICOM files stored locally on the user's computer (Figure 2).
   NOTE: Inputs to the module M-mode view DICOM image analysis from a specific patient and LV parameter extraction are the testing image (provided by the user in DICOM format representing M-mode or Apical view) and static files (Figure 1). The static files comprise a) the input.txt file, which is a template file used to overwrite the calculated output measures from the module (new output files are created, in turn, from this file), and b) the template image file, which is used for template matching (i.e., template_dicom_GEMS.jpg).
4. Select either Private or Public folders as the Destination Folder for the files. Type the desired comment or note in the Comments section before starting the workflow. Select the Apical view and M-mode LV ultrasound images and DICOM files uploaded previously (Step 1.3.) (Figure 3).
   NOTE: Ultrasound Heart Segmentation tool consists of two sub-modules: the apical view and the M-mode view. Both modules are accessed via the platform. The analysis includes the calculation of the characteristic parameters visible in corresponding cross-sections. In both cases, the necessary values are written to a file that is taken as input to the parametric heart model from the user side (see the discussion for details).
5. Click on the Execute button. The interactive platform notifies the user when the running workflow is finished.
6. Visualize the created geometry of LV directly on the platform (a 3D model will appear automatically on the screen, which can be rotated using the mouse). Available options include Shaded and Wireframe models for visualization.
   ​NOTE: After the LV dimensions are extracted, based on mesh options, the mesh of finite elements is generated, and finite element analysis simulation is run with user-prescribed inlet and outlet velocities (see below).

2. PAK finite element solver for fluid-structure (FS) simulations

NOTE: This tool can be used for the finite element analysis of the coupled solid-fluid problems. It supports both strong and loose coupling between solid and fluid. Elements can be hexahedrons or tetrahedrons, with or without an additional node at the center of the element. This solver has built-in material models such as the Holzapfel model, hunter muscle model, etc. The PAK-FS information flow diagram is depicted in Figure 4. It starts from the input file and PAK preprocessor. The PAK preprocessor tool outputs a DAT file, which will serve as an input file for the finite element solver. The final output from the solver is VTK files that include the results of the finite element simulation: velocities, pressures, strains, and stresses in the left ventricle.

1. Download the template files for the boundary conditions of inlet and outlet velocities using the buttons in the bottom section (Figure 3); if patient-specific flow boundary conditions are available, download and use these files. Download the mesh options by clicking on the corresponding button in the bottom section. Save these files into either Private or Public folders.
2. Upload these files in a similar way as uploading images (Figure 3). The prescribed inlet and outlet velocities simulate the drug condition, while the mesh options control the density of the finite element mesh. For simulating patient-specific conditions, modify the default values of pressure, flow, material properties, and calcium function.
3. Click on the Execute button. A new running workflow will appear in the list. If any of the sections of the workflow are not clear, click on the Help File button (Figure 3) to view detailed instructions on how to use this workflow, as well as to interpret the results.
   NOTE: If everything executes without errors, the status of the workflow will change from "running" to "finished ok".
4. View the results with several options. Alternatively, download the results; the results folder contains VTK files, CSV files, and animations.
   1. Click on the eye button to view ejection fraction and global work efficiency values, as well as diagrams of pressure vs. volume, pressure vs. strain, and myocardial work vs. time.
   2. Click on the camera button to preview and play animations of the displacement, pressure, shear stresses, and velocity fields.
   3. Click on the 3D visualization button to visualize the outputs online in ParaView Glance.
      1. Load multiple VTK files previously downloaded as results. See several parameters of interest and change the field, for instance, to velocity for visualization.
      2. Rotate the model or change the color scheme. Choose Surface with edges or Wireframe for the representation of the surface. Apply the same methodology to every loaded VTK file.

3. Determination of ventricular activation sequence from ECG measurement

NOTE: A modified FitzHugh-Nagumo model of the cardiac cell was implemented. The precordial leads were modeled with the standard six electrodes. The potential of the heart was optimized with inverse ECG. Starting from activation in the sinoatrial node (which is a function of time), with heterogeneous action potential through the heart and torso, the user can get electrical activity on the total torso model. The PAK-TORSO information flow diagram is given in Figure 5. The user provides scaling of the torso model in all directions (x,y,z) and ECG signal function. Then the scaled model is created, and its behavior is simulated using the PAK-FS solver. The user provides these input values in a text file. The output of the simulation is a VTK file with the electrical activity of the heart in the torso-embedded environment.

1. On the homepage, go to Execute Workflow and then choose torse-cwl in the list of available workflows. Add a comment or note in the Comments section before executing the workflow.
2. Click on the Input Template File button and save the content shown on the webpage as input.txt file, which will be used for the torso model. In the Input File field, select the downloaded input.txt file. After the file is imported, click on the Execute button to start the calculation.
   NOTE: A new workflow will appear in the bottom left, with its status shown as "running"; the calculation lasts about 0.5 h, and then the status will change to "finished ok".
3. Click on the eye or camera buttons in the bottom left corner to visualize the available simulation reports or animations directly on the platform.
   NOTE: The results include distributions of the electrical field, velocity, pressure, shear stress, and deformation at each time step. Animations for each of these distributions are possible to view for the total heart cycle.
4. Alternatively, click on the 3D visualization button to visualize the outputs online in ParaView Glance.
   1. Select the OPEN A FILE button, go to the GIRDER tab, enter user credentials if prompted, and open the Private folder.
   2. On the next page, select the workflow-outputs folder, and open the torso-cwl folder. Open the first folder on the list.
   3. See the list of VTK files that represent the simulation results. Choose one or more files and click on the SELECT button to load the file in ParaView Glance.
   4. Manipulate the model geometry (i.e., move, rotate, zoom in or out, etc.) using the mouse.
   5. Select different options for model representation as follows.
      1. Choose the Wireframe option to see the interior of the torso with a heart incorporated within the torso. Choose the Points option to display a dotted representation of the torso model with full heart mesh.
      2. Adjust the Point Size value to change the display results. Adjust the Opacity value to see the interior of the torso and display results inside the heart mesh.
      3. Click on the Color By drop-down menu and choose the desired option, for instance, electric potential. Change the default color scale to any of the listed options.

As an example, the workflow for ultrasound M-mode DICOM image analysis and LV parameter extraction is presented in Figure 1. M-mode and apical view could be tested separately or one after another, depending on the parameters of interest. If tested one after another, the results are appended to one common file (separately for systole and diastole phases). If only one view is tested, the values for unknown parameters are taken from the default file input.txt (Figure 1). The SILICOFCM tool analyzes the images depending on the view mode, as follows.

Apical view analysis includes segmenting the LV using a previously trained U-net framework and calculating the bordering rectangle, based on which the parameters LV length, diastolic, 2D (LVLd[cm]) and LV length, systolic, 2D (LVLs[cm]) are calculated. The user should define if the view represents the systolic or diastolic phase. The output files are as follows: a) input_parametric_diastole.stl or input_parametric_systole.stl, which represents the input file for PAK solver with created geometry, including the calculated value for LVLd[cm] or LVLs[cm] depending on the user input of whether the image represents diastole or systole; b) output_Apical_view.txt, which contains the data for LVLd[cm] or LVLs[cm], depending on the user input of whether the image represents diastole or systole.

M-mode view analysis includes the bordering of the characteristic areas of LV in this view. Based on these areas, various parameters are calculated. The output files are as follows: a) input_parametric_diastole.stl and input_parametric_systole.stl, which represent the input geometry for PAK solver, including calculated values for interventricular septum thickness, diastolic, M-mode (IVSd[cm]), interventricular septum thickness, systolic, M-mode (IVSs[cm]), LV internal dimension, diastolic, M-mode (LVIDd[cm]), LV internal dimension, systolic, M-mode (LVIDs[cm]), LV posterior wall thickness, diastolic, M-mode (LVPWd[cm]), and LV posterior wall thickness, systolic, M-mode (LVPWs[cm]); b) output_M_mode.txt, which contains the data for IVSd [cm], IVSs [cm], LVIDd [cm], LVIDs [cm], LVPWd [cm], and LVPWs [cm].

To generate an example of LV, as realistic as possible, a parametric type of model with specific parts, namely the base part, valves (aortic and mitral), and connecting part (i.e., the connection between base and valves), was generated via PAK finite element simulations (Step 2) using the calculated parameters of LV geometry. Every part has its length and radius, as well as the number of layers (division). The finite element model consists of a fluid domain surrounded by a solid wall that is connected to the base and part of the connection layer. To simulate the realistic behavior of the LV model, prescribed functions for velocities at the inlet (mitral; Figure 6A) and outlet (aortic; Figure 6B) valves were used. An algorithm for the automatic calculation of fiber direction was applied to this finite element model. The results for one layer and three-layered solid representation are shown in Figure 7.

Figure 8 displays the pressure distribution inside the parametric model during a time cycle of 1 s, divided into 10 time steps. During the first five steps, there are no significant pressure changes until the contraction, when it reaches the highest value on the scale and the volume of the model starts to increase at the same time. Figure 9 shows how the velocity is distributed inside the fluid part of the LV model. As shown, there are notable value peaks inside the branches, caused by fluid flow during the loading/unloading cycle. Figure 10 shows how the displacements are distributed along the model. Similar to the pressure change, during the first two steps, the displacements are negligible until the contraction, when they become maximal at the bottom part of the model. Throughout the remaining time, the model slowly returns to its undeformed state. Figure 11 shows the resulting pressure-volume (PV) diagram. Figure 12 shows the vectorial representation of velocities inside the fluid part of the ventricle model. Phase velocity vectors are present inside the mitral valve during the fluid inflow and, after the ventricle is filled with fluid from 0.7 s, fluid goes out through the aortic part.

Monodomain equations for the fully coupled heart torso system were applied in order to simulate electrophysiology. Myocardial and torso conductivities are defined in Table 1^{10,11,12,13,14}. Zero flux was used for V_m as a boundary condition for all the interior boundaries in contact with the lungs, torso, and cardiac. Therefore, −n · Γ = 0, where n is the unit outward normal vector on the boundary and Γ is the flux vector through that boundary for the intracellular voltage, equal to Γ = − σ· ∂V_m/∂_n. Next, classical approaches were implemented for solving the ECG inverse problem using the epicardial potential formulation. The methods used were the family of Tikhonov methods and L regularization-based methods^{10,11,12,13,14}.

Whole-heart electrical activity in the torso-embedded environment is presented in Figure 13. A specialized conduction system with heterogeneous action potential morphologies throughout the heart was incorporated in the initiation of activation in the sinoatrial node. Figure 14 shows the maps of body surface potential in a healthy subject during the progression of ventricular activation in nine sequences that correspond to the measured ECG signal.

Figure 1: Information flow diagram for Ultrasound M-mode DICOM image analysis and LV parameter extraction. Please click here to view a larger version of this figure.

Figure 2: Uploading new files for analysis. Please click here to view a larger version of this figure.

Figure 3: User interface for Ultrasound Apical view DICOM image analysis and LV parameter extraction. Please click here to view a larger version of this figure.

Figure 4: PAK-FS information flow diagram. Please click here to view a larger version of this figure.

Figure 5: PAK-TORSO information flow diagram. Please click here to view a larger version of this figure.

Figure 6: Velocity profile. (A) Inlet function of velocity at mitral valve cross-section, and (B) outlet velocity function at aortic valve cross-section. Please click here to view a larger version of this figure.

Figure 7: Parametric LV model with fibers direction. (A) One-layered solid wall, and (B) three-layered solid wall representation. Please click here to view a larger version of this figure.

Figure 8: Pressure field inside the parametrical LV model. Five different time steps are presented. Please click here to view a larger version of this figure.

Figure 9: Velocity field inside the parametrical LV model. Five different time steps are presented. Please click here to view a larger version of this figure.

Figure 10: Displacements in the parametrical LV model. Five different time steps are presented. Please click here to view a larger version of this figure.

Figure 11: Pressure vs. volume (PV) diagram for LV fluid-structure interaction model. The prescribed cycle duration of 1 s is divided into 10 time steps. Five representative time steps are shown. Please click here to view a larger version of this figure.

Figure 12: Vectorial representation of velocities in the parametrical model of the left ventricle. Four representative time steps are shown. Please click here to view a larger version of this figure.

Figure 13: Simulation of whole-heart activation at various time points on the lead II ECG signal. 1–9 activation sequences in (B) correspond to the ECG signal in (A). The transmembrane potential in mV is denoted by the color bar. Please click here to view a larger version of this figure.

Figure 14: Maps of body surface potential in a healthy subject. Progression of ventricular activation in nine sequences (lower panel) corresponding to the ECG signal (upper panel). The heart activity range in mV is denoted by the color bar. Please click here to view a larger version of this figure.

Table 1: Parameters for monodomain model with modified FitzHugh-Nahumo equations.

The SILICOFCM project is an in silico clinical trials platform to design virtual patient populations for risk prediction, testing the effects of pharmacological treatment, and reducing animal experiments and human clinical trials. Testing the effects of pharmacological treatment was modeled with prescribed inlet/outlet boundary flow conditions, calcium function, and material wall properties. This platform integrates multiscale methods on the sarcomeric level with whole-heart performance and the functional optimization level for the risk prediction of patient-specific conditions during cardiomyopathy disease progression.

Heart modeling for cardiomyopathy and electromechanical coupling of the LV in the SILICOFCM⁸ project has been presented. The geometry of the heart model included seven different regions: 1) sinoatrial node; 2) atria; 3) atrioventricular node; 4) His bundle; 5) bundle fibers; 6) Purkinje fibers; and 7) ventricular myocardium. Body surface potential maps in a healthy subject during the progression of ventricular activation in nine sequences corresponding to the ECG signal have been presented.

After the dimensions of the left ventricle are extracted, based on mesh options, the mesh of finite elements is generated and finite element analysis simulation is run with user-prescribed inlet and outlet velocities. Users can visualize solutions directly on the platform by looking at the available animations and diagrams. Users can visualize pressure-volume, pressure-strain, and myocardial work-time diagrams by clicking on the eye symbol. If users click on the camera button, the list of available animations of various fields (displacements, pressures, velocity, shear stresses) appears.

There are some limitations to the study. Parametric LV geometry was extracted from ultrasound images. The future version will go into more detailed geometrical reconstruction. A direct connection with specific drugs for cardiomyopathy is not presented in this manuscript. In this study, it is controlled through boundary conditions for flow and pressure. Genetic data are currently incorporated through calcium heart muscle function and the nonlinear material property of the wall. In the future version, more details from image reconstruction and genetic data for patient-specific cardiomyopathy disease (hypertrophic and dilated genetic variant) will be considered.

The SILICOFCM computational platform⁸ will open a new avenue for in silico clinical trials, specifically for cardiac disease and risk prediction for the patient-specific condition. The gold standard for today's clinical practice for risk prediction is the standard survival calculator for cardiomyopathy patients. The platform used here can give more information in comparison with the current medical standard, as the modeling includes not just the biomarkers but also patient-specific geometry, flow and pressure hemodynamics conditions, wall material properties (from deformation in the images), and drug response with different combinations of boundary conditions, calcium function, and material properties.