PATIENT-SPECIFIC PREDICTION OF ABDOMINAL AORTIC ANEURYSM EXPANSION USING EFFICIENT PHYSICS-BASED MACHINE LEARNING APPROACHES By Zhenxiang Jiang A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Engineering Mechanics—Doctor of Philosophy 2022 ABSTRACT PATIENT-SPECIFIC PREDICTION OF ABDOMINAL AORTIC ANEURYSM EXPANSION USING EFFICIENT PHYSICS-BASED MACHINE LEARNING APPROACHES By Zhenxiang Jiang Computational vascular Growth and Remodeling (G&R) models have been developed to capture key physiological and morphological features during the arterial disease progression and have shown promise for aiding clinical diagnosis, prognosis prediction, and staging clas- sification. However, the translation of computational G&R models into their applications has yet to wait for clinical practice. Partly, due to the high complexity of the arterial adap- tation mechanism, high-fidelity arterial G&R simulations typically require hours or even days, which hinders its time-consuming applications such as patient-specific parameter es- timation, disease prediction, verification, validation, and sensitivity analysis. Furthermore, the typical Finite Element Method (FEM) based computational G&R model should be ex- tended to provide the uncertainty quantification associated with simulation and prediction results. Therefore, to enhance practicality of the G&R modeling, we develop a novel and computationally efficient simulation framework that comprehensively combines physics-based G&R simulations and data-driven machine learning methods using a Multi-Fidelity Surro- gate (MFS) approach. This greatly enhances the computational efficiency of arterial G&R simulations, enabling more time-consuming applications such as personalized parameter es- timation. The proposed framework is then tested for a specific disease, Abdominal Aortic Aneurysms (AAAs), by estimating G&R model parameters from follow-up CT images in 21 patients. The physic-based machine learning simulation framework of G&R enables computationally efficient personalized simulations that underlie the prediction of patient-specific progression of AAAs within clinically relevant time frame (1∼2 hours). To the end, two different physics- based machine learning techniques are developed to predict the evolution of AAA geometries. First, by employing a cokriging method, a set of follow-up CT images from AAA patients and the corresponding personalized G&R simulations are integrated to predict their patient- specific evolutions. Second, G&R simulations are augmented into a massive in silico dataset to capture variations in real patient AAA geometries. This massive in silico dataset, along with a limited medical follow-up dataset, is used to train a Deep Belief Network (DBN) to provide fast predictions of patient-specific AAA expansion. These efforts not only demon- strate the vital role of physics-based machine learning methods and computational G&R models in predictive medicine, but also provide insights into the creation of physics-based digital twins of humans. Accordingly, at the end of this dissertation, the proposed frame- work is further extended to generate Virtual Patient Cohort (VPC) of a targeted disease population, which is illustrated for a simplified but representative model of mitral valve regurgitation. Copyright by ZHENXIANG JIANG 2022 ACKNOWLEDGMENTS First of all, I would like to express my sincere thanks to my advisor, Dr. Seungik Baek, for his dedicated support and professional guidance. This research would not have been possible without his immense knowledge, insightful feedback and patience. I would also like to express my gratitude to Michigan State University and National Institutes of Health for their generous sponsorship of my research. Besides my advisor, I thank the other members of my dissertation committee, Dr. Lik Chuan Lee, Dr. Zhaojian Li and Dr. Yuying Xie for their insightful comments and advises. I also thank my colleagues in the Cardiovascular and Tissue Mechanics Research Laboratory for enlightening discussions and good days over the past years. My sincere thanks also goes to Tom Battisti, Dr. Jiang Yao and Dr. Steve Levine, for offering me the great internship opportunity for the ENRICHMENT in Silico Clinical Trial Project sponsored by Dassault Systèmes and U.S. Food and Drug Administration. Last, I would like to thank my girlfriend Yimin Wu for her warm companion over the past four years, especially during the pandemic. I am grateful to my parents, Qikun Jiang and Wanhong Liu, as well as my grandmother, Shuying Sheng, who brought me to this wonderful world and supported me spiritually through my life. v TABLE OF CONTENTS LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x LIST OF ALGORITHMS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii Chapter 1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Arterial growth and remodeling models and limitations . . . . . . . . . . . . 2 1.3 Data-driven predictions for arterial diseases . . . . . . . . . . . . . . . . . . 4 1.4 Objective . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Chapter 2 AN EFFICIENT REPRODUCTION OF STRESS-MEDIATED VASCU- LAR GROWTH AND REMODELING MODELS . . . . . . . . . . . . . . . 7 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2 Mathematical modeling of vascular G&R . . . . . . . . . . . . . . . . . . . . 9 2.2.1 Kinematics and configurations . . . . . . . . . . . . . . . . . . . . . . 9 2.2.2 Constitutive laws and constituent turnover . . . . . . . . . . . . . . 12 2.2.3 Homeostasis and stress-mediated adaption . . . . . . . . . . . . . . . 14 2.2.4 Elastin degradation and smooth muscle active tone . . . . . . . . . . 15 2.3 High-fidelity 3D membrane G&R model . . . . . . . . . . . . . . . . . . . . . 17 2.3.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.3.2 Membrane element and local coordinate system . . . . . . . . . . . . 19 2.3.3 Finite element formulations . . . . . . . . . . . . . . . . . . . . . . . 20 2.3.4 Code verification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.4 Low-fidelity axisymmetric G&R model and verification . . . . . . . . . . . . 24 Chapter 3 PHYSICS-BASED SURROGATE MODELS . . . . . . . . . . . . . . . . . . . 27 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.2 Probabilistic collocation method . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.2.1 Deterministic input . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.2.2 Stochastic input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.2.3 Selection of base functions and collocation points . . . . . . . . . . . 31 3.3 Multi-fidelity surrogate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.3.1 Kriging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.3.2 Cokriging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.4 Demonstrations of multi-fidelity surrogate . . . . . . . . . . . . . . . . . . . 40 3.4.1 A univariate demonstration of MFS . . . . . . . . . . . . . . . . . . . 41 vi 3.4.2 A MFS-based approximation to the G&R models of AAA growth . . 42 3.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 Chapter 4 AN AUTOMATIC PATIENT-SPECIFIC PARAMETER ESTIMATION FRAMEWORK USING MULTI-FIDELITY SURROGATE . . . . . . . . . 46 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 4.2 Patient-specific parameter estimation framework . . . . . . . . . . . . . . . . 50 4.2.1 Data assimilation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.2.2 Objective function and uncertainty quantification . . . . . . . . . . . 52 4.2.3 Lower confidence bound-based optimal sampling . . . . . . . . . . . . 53 4.3 Case study setup and clinical data . . . . . . . . . . . . . . . . . . . . . . . . 54 4.3.1 Clinical data of AAAs . . . . . . . . . . . . . . . . . . . . . . . . . . 54 4.3.2 G&R models and MFS setup . . . . . . . . . . . . . . . . . . . . . . 56 4.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 Chapter 5 A PHYSICS-BASED DEEP LEARNING APPROACH TO PREDICT AB- DOMINAL AORTIC ANEURYSM ENLARGEMENT . . . . . . . . . . . . . 65 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 5.2 Deep neural network algorithm . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.3 Data processing and results . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 5.3.1 Data processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 5.3.2 Test set-up . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 5.4 Test prediction results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.5 Monte-Carlo cross-validation . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 Chapter 6 A PHYSICS-BASED MACHINE LEARNING APPROACH TO PREDICT ABDOMINAL AORTIC ANEURYSM ENLARGEMENT . . . . . . . . . . 82 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 6.2 Personalized prediction of AAA extension . . . . . . . . . . . . . . . . . . . 83 6.3 Results of AAA expansion prediction . . . . . . . . . . . . . . . . . . . . . . 85 Chapter 7 A VIRTUAL PATIENT ENGINE FOR THE CREATION OF A MITRAL VALVE REGURGITATION PATIENT COHORT . . . . . . . . . . . . . . . 88 7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 7.2 Material parameter calibration . . . . . . . . . . . . . . . . . . . . . . . . . . 90 7.3 Parameterization of mitral valve apparatus geometry . . . . . . . . . . . . . 91 7.4 Automatic simulation framework . . . . . . . . . . . . . . . . . . . . . . . . 93 7.5 Design of experiment setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 7.6 Virtual patient engine . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 vii 7.6.1 Create surrogate model . . . . . . . . . . . . . . . . . . . . . . . . . . 97 7.6.2 Generate surrogate-based VPC . . . . . . . . . . . . . . . . . . . . . 98 7.6.3 Down select VPC for a targeted population . . . . . . . . . . . . . . 99 7.6.4 Find representative VPs using a clustering algorithm . . . . . . . . . 101 7.7 Results: physics-based VPC and statistical analysis . . . . . . . . . . . . . . 102 7.8 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 viii LIST OF TABLES Table 3.1 Comparison between two-fidelity G&R models. . . . . . . . . . . . . . . . 41 Table 4.1 Demographic data and CT scan information of individual patients. . . . . 55 Table 4.2 Information of patient-specific parameter optimization. . . . . . . . . . . . 60 Table 5.1 Demographic data of patients . . . . . . . . . . . . . . . . . . . . . . . . . 73 Table 5.2 Effect of the number of nodes in a 2-layer DBN on the model testing. . . 76 Table 5.3 The absolute and relative prediction errors of 6 testing samples under DBN and Mix-effects model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 Table 6.1 Prediction results of personalized AAA expansion. . . . . . . . . . . . . . 86 Table 7.1 Results of MVA material parameter estimation. . . . . . . . . . . . . . . . 91 Table 7.2 List of input parameters and their ranges. . . . . . . . . . . . . . . . . . . 95 Table 7.3 A demonstration of LHS samples with the values of input parameters and targeted outputs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 Table 7.4 Performance of different surrogate approaches in predicting anatomic orifice area (mm2 ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 ix LIST OF FIGURES Figure 2.1 Important configurations of vascular G&R, including the current in-vivo configurations κτ , the computational reference configuration κR , and the natural stress-free configuration of the newly produced constituent κkn (τ ), in which time τ ∈ [0, t]. The deformation gradient F(τ ) maps from κR to κτ . Gk (τ ) denotes the deposition pre-stretch mapping from κkn (τ ) to κτ . . . . . . . . . . . . . . . . 10 Figure 2.2 X = {X1 , X2 , X3 } are global Cartesian coordinates defined on base vec- tor {E1 , E2 , E3 } in reference configuration. {Ee1 , Ee2 , Ne } are transformed base vectors that define the local coordinates. The centroid Xc and nodal coordi- nates of elements {X(1) , X(2) , X(3) } are defined by global coordinates. The asso- ciated nodal coordinates in current configuration of G&R process are provided by {x(1) , x(2) , x(3) }. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 Figure 2.3 A verification between two 3D membrane G&R simulations implemented on FEniCS and MATLAB, respectively. The same simulation setup and input parameters are used in two simulations. . . . . . . . . . . . . . . . . . . . . . . . 24 Figure 2.4 Two-fidelity G&R FEM models. The 3D G&R model is built on a three dimensional geometry of artery with linear membrane elements. The axisymmet- ric G&R model is built on an idealized axisymmetric geometry with quadratic membrane elements. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 Figure 2.5 Verification results of AAA growth between the FEniCS-based high-fidelity model (represented by ‘FEniCS’) and the MATLAB-based low-fidelity model (represented by ‘Matlab’). The maximum diameter (max. d) vs. time (t) for both simulations are plotted in Figure a. The diameter curves (d) along the length of arterial centerline (s) for both simulations at t = 2000 days and t = 3000 days are plotted in Figure b. Note that the crosses and circles in Figure b represent the AAA diameters at nodal locations. . . . . . . . . . . . . . . . . . . 26 Figure 3.1 A simple demonstration to cokriging-based MFS. s and d denote the length along the centerline and the associated diameter, respectively. The red and blue dashed lines represent the physics-based simulation outputs from HFM and LFM, respectively. The kriging-based surrogate is trained from 6 HF samples while the cokriging-based MFS is trained from the same set of 6 HF samples as well as 50 additional LF samples. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 Figure 3.2 Contour plots of the approximated maximum diameter and the associated uncertainty over the map of Kg and σdmg at t = 4000 days. The red dots represent the HF samples. The unit of color bar is mm. . . . . . . . . . . . . . . . . . . . 43 Figure 3.3 The MFS approximations of diameters against independent input vari- ables in which irrelevant input parameters are fixed. The approximated error is indicated in blue shades. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 x Figure 4.1 Schematic flow of the proposed MFS-based parameter estimation algorithm. 49 Figure 4.2 A contour plot that provides the LCB surface Clcb (θ) at the final iteration of optimization. The unit of color bar is mm. The red dots denote the additional samples used to iteratively update the MFS. The black arrows indicate the order of additional samples. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 Figure 4.3 Diameter (d) versus length along centerline (s) of the AAA of Patient 18. The red dashed lines with circles represent the clinical data measured from follow-up CT images captured at t = [0, 1.68, 2.87, 4.15] years. The black curves represent MFS-based patient-specific G&R simulation outputs. The uncertainty related to MFS approximation is indicated by 95% confidence interval and is shown by blue shades. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 Figure 5.1 Diagram of overall methodology. A massive number of in silico data and a small number of patient data are collected, which is followed by a two-stage model training and a model testing. . . . . . . . . . . . . . . . . . . . . . . . . . 68 Figure 5.2 The deep architecture of the DBN. Two layers of RBM are trained in an unsupervised manner (pre-trained) using CD-1 algorithm. The top layer utilizes a neural network sigmoid regression for the prediction. . . . . . . . . . . . . . . 69 Figure 5.3 Diagram of the model training. The DBN is pre-trained by in silico data in an unsupervised manner, and is unfolded into a neural network. Next, the neural network is fine-tuned with the patient data in a supervised manner to predict the AAA expansion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 Figure 5.4 Effect of the number of epochs on the prediction error. The RMSE and the fine-tuning training time are plotted in solid blue and dashed red lines, re- spectively. The training time increases linearly with the number of epochs, while the RMSE rapidly decreases at the beginning but converges at around 400 epochs. 77 Figure 5.5 The true value, the DBN model prediction and mixed-effects model pre- diction of IMDCs are shown in dotted dashed black (‘true’), solid red(‘DL pre- diction’), and dashed green lines(‘ME prediction’), respectively. s denotes the coordinate of location along the centerline of AAA and d represents the associ- ated maximum diameter measured by inscribed sphere method. . . . . . . . . . 79 Figure 6.1 Overall schematic flow of the proposed framework. A MFS is approxi- mated from G&R models using a cokriging method. The in silico data of patient- specific simulations are obtained from both the MFS and the clinical data using a parameter estimation algorithm. The clinical and the in silico data are combined to predict the progression of arterial disease. . . . . . . . . . . . . . . . . . . . . 84 xi Figure 6.2 Two examples that compare prediction results with the ground truth of AAA shape curve, where diameter (d) is plotted against the length along cen- ter line (s). The blue shades represent the 95% confidence interval of AAA growth prediction. The black circles and red crosses represent the ground truth of aneurysmal diameters that locate within and outside the confidence range, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 Figure 7.1 A demo simulation of parameterized mitral valve apparatus model. . . . 93 Figure 7.2 A diagram of the automatic simulation framework for parameterized MVA model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 Figure 7.3 Simulation outputs vs. surrogate outputs of different surrogate models. . 97 Figure 7.4 Simulation outputs vs. surrogate output of other target quantities, in- cluding tenting angle of the anterior leaflet, the tenting angle of posterior leaflet and the tenting height. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 Figure 7.5 Histogram of orifice area for surrogate-based VPs. . . . . . . . . . . . . . 99 Figure 7.6 Histogram of orifice area for surrogate-based VPs and the down selection criterion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 Figure 7.7 A demonstration figure with 200 Monte-Carlo samples displayed on a two- dimensional parameter space of Bf vs. leaflet ratio (LLR). Included samples are indicated by red dots, and the excluded samples are indicated by blue circles. . . 101 Figure 7.8 Included Monte-Carlo samples and 50 VP samples in a 3D parameter space.103 Figure 7.9 Top-down view and A2P2 cut view of 12 VPs in the VPC. . . . . . . . . 104 Figure 7.10 A histogram of the orifice area extracted from 50 virtual patients in the physics-based VPC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 xii LIST OF ALGORITHMS Algorithm 1 The Probabilistic Collocation Method (PCM) model. . . . . . . . . 34 Algorithm 2 MFS-based parameter estimation for arterial G&R models . . . . . . 50 xiii Chapter 1 INTRODUCTION 1.1 Motivation Arterial diseases, such as Abdominal Aortic Aneurysm (AAA) and thoracic aortic aneurysm, are major risks to human health and pose a substantial social and economic burdens. These diseases vary by arterial types and pathological mechanisms, but share similarities in irregu- lar turnover of constituents and changes of mechanical and geometrical properties. Therefore, capturing these pathological and morphological changes in advance is crucial for early diag- nose and control of arterial disease progression, which requires both good understandings of arterial diseases and sophisticated computational techniques. To achieve this, we developed computational vascular Growth and Remodeling (G&R) mod- els [Baek et al., 2006, Zeinali-Davarani et al., 2011], which enable patient-specific simulation of arterial diseases progression by accounting for changes and adaptions in arterial tissue. However, the translation of computational G&R models into their applications has yet to wait for clinical practice. This is mainly because the high complexity of the arterial adap- tation mechanism limits the computational efficiency of the computational G&R models 1 and hinders its time-consuming applications such as patient-specific parameter estimation, disease prediction, sensitivity analysis and validation. Therefore, in this study, we aim to develop a computationally efficient simulation frameworks that combine physics-based G&R simulations and data-driven machine learning methods for time-consuming clinical applica- tions. 1.2 Arterial growth and remodeling models and limi- tations Cardiovascular models fueled by multifidelity data and fundamental laws of biomechanics have been widely applied for clinical tasks [Borghi et al., 2008, Szafron et al., 2019, Taylor et al., 1999]. Among these models, the stress-mediated arterial G&R model [Baek et al., 2006, Cyron et al., 2016, Humphrey and Rajagopal, 2002] has been proposed as an essential tool for capturing important features of anatomical, biological, and pathological changes, as well as predicting disease progress, thus serving clinical applications such as diagnosis, stratification, and disease treatment. Its potential has been demonstrated under various pathological situations, such as AAAs [Zeinali-Davarani et al., 2011], thoracic aortic aneurysms [Mousavi et al., 2019], cerebral aneurysms [Baek et al., 2005], stent implantation [Laubrie et al., 2019], and pulmonary autografts [Famaey et al., 2018]. In this study, we test our new framework primarily for a specific arterial disease, abdominal aortic aneurysm (AAA). An AAA is a localized enlargement of the abdominal aorta that undergoes gradual dilation over years and can eventually lead to a fatal rupture[Klink et al., 2011, Vardulaki et al., 1999]. Typically, abdominal aortas with maximum diameters larger than 3 cm are considered as AAA, and the growth rate of AAA diameter ranges from 0.1 cm/year to 1.5 cm/year depending on the cohort groups [Bhak et al., 2015, Matsuda et al., 2 1992, Vega de Céniga et al., 2006]. Due to the high mortality (80%) of rupture [Thompson et al., 2012] , effective factors for estimating rupture risk have been collected, including growth rate [Brown et al., 2003, Limet et al., 1991], thrombus volume [Parr et al., 2010], aortic volume [Martufi et al., 2013], asymmetry [Scotti et al., 2005] and tortuosity [Cappeller et al., 1998]. In addition, recent studies (Akkoyun et al. [2020], Lee et al. [2017b], Siika et al. [2015]) suggested that the growth rate of AAA is the key predictor of AAA enlargements. These findings motivated us to develop computational models which can capture the mechanical properties of AAA and thus provide better predictions of longitudinal growth of AAA. The arterial G&R model is able to capture the geometrical changes of AAAs in response to external stimuli and changes of environment, and thus simulate the localized aortic wall bulge, which is the most significant feature of AAAs. In addition, other geometrical and biological features of AAAs such as asymmetry, tortuosity, accumulated thrombus, and lost of elastin and smooth muscle, have been proved to play essential roles in AAA extension [Akkoyun et al., 2020, He and Roach, 1994, Kwon et al., 2015]. For example, the reduction of elastin and smooth muscle causes a compensatory mechanism by which a stress-mediated collagen accumulation protects the aortic wall, but, on the other hand, the accumulated collagen reduces the extensibility of aortic wall. As such a way, the comprehension of elastin degradation, collagen turnover and their biomechanical properties such as production rate, half-life, deposition stretch and stiffness are associated with the AAA growth and wall weak- ening in predicting risk of rupture [Wilson et al., 2013]. These features of AAAs have been successfully captured in recent studies of stress-mediated G&R models, and have been successfully demonstrated in feasibility studies [Baek et al., 2006, Volokh and Vorp, 2008, Watton and Hill, 2008, Zeinali-Davarani et al., 2011]. However, the potential of clinical oriented applications for such G&R models is still hindered by various factors such as (1) the deterministic nature of physics-based G&R models that restricts the use of uncertainty analysis, (2) the difficulties of taking inter-personal patholog- 3 ical diversity into account, and (3) the low computational efficiency attributed by the high complexity of arterial diseases. Especially, the third factor, i.e. low computational efficiency, hampers the implementation of time-consuming applications of vascular G&R models such as patient-specific parameter estimation, disease prediction, verification, validation, and sen- sitivity analysis. Especially, Ambrosi et al. [2019] suggested that one of the next priority in the biomedical engineering community should be significantly improving the computational efficiency to enable personalized modeling. This can be done through reducing the high complexity of fundamental laws of G&R and using methods such as data-driven machine learning. 1.3 Data-driven predictions for arterial diseases Data-driven approaches have been explored to enhance the prediction capability on the progress of vascular diseases and their risk. Take aortic aneurysm as an example. Jordanski et al. [2018] implemented multiple machine learning models to predict the wall shear stress of AAAs . Hirata et al. [2020] utilized an eXtreme Gradient Boosting (XGBoost) approach to predict the growth and classify fast growth of AAAs by selecting different features. Do et al. [2019] developed the Dynamical Gaussian Process Implicit Surface (DGPIS) modeling to predict AAA growth, using longitudinal computer tomography (CT) scans of AAAs that are captured at different times in a patient-specific way. Zhang et al. [2019] developed a G&R model based Bayesian calibration approach that predicted the AAA enlargement of three patients. Nevertheless, these models lack either sufficient computational efficiency or the implementa- tion of fundamental biomechanical laws. In this study, we therefore aim to develop an inno- vative strategy for a computational arterial G&R modeling framework that integrates both physics-based simulations and data-driven machine learning approaches, i.e. a physics-based 4 machine learning framework. This novel framework is expected to enhance the computational efficiency while taking into account of both uncertainty and inter-personal variability, thus providing effective and efficient predictions of arterial diseases within the clinical relevant computational time frame (1-2 hours) [Taylor et al., 1999]. Furthermore, the proposed physics-based machine learning framework holds the potential to generate digital twins of Virtual Patient Cohort (VPC), which generally allow the physics- based simulations to capture the inter-patient variety in a targeted cardiovascular disease [Niederer et al., 2020]. This mainly benefits from the high computational efficiency of the proposed physics-based machine learning framework. In this study, therefore, we extend the proposed framework to allow a fast workflow that can rapidly generate a large number of approximate simulations and down-select them into a small number of models to represent VPC. This workflow is considered a Virtual Patient Engine (VPE), and is illustrated for a simplified but representative model of mitral valve regurgitation. 1.4 Objective In general, this study aims to develop an efficient physics-based machine learning framework that comprehensively combines G&R models and machine learning approaches to enable time-consuming clincal applications. Toward this goal, specific aims are set as: • Implement G&R theories on FEniCS to reproduce an efficient simulation framework of AAA growth (Chapter 2). • Create surrogate models to approximate time-consuming G&R simulations using ma- chine learning (Chapter 3). • Develop a computationally efficient framework to estimate patient-specific parameters of AAAs using multi-fidelity surrogate (Chapter 4). 5 • Develop a deep learning predictive tool to predict AAA growth using massive G&R in silico data and a limited clinical longitudinal data (Chapter 5). • Develop a physics-based machine learning approach to predict AAA enlargement from patient-specific simulations and clinical data (Chapter 6). • Extend the proposed physics-based machine learning framework to create a workflow (virtual patient engine) that generates virtual patient cohort for a target disease popu- lation (Chapter 7). This workflow is demonstrated by a cardiac disease, i.e., secondary mitral regurgitation. 6 Chapter 2 AN EFFICIENT REPRODUCTION OF STRESS-MEDIATED VASCULAR GROWTH AND REMODELING MODELS 2.1 Introduction The vascular G&R theory has been employed to model the stress-mediated adaptation of vascular tissues to capture long-term changes of geometrical features in AAAs [Baek et al., 2006, Volokh and Vorp, 2008, Watton and Hill, 2008]. More specifically, in our G&R models, the aortic tissue is considered to be a pre-streched object that adapts to inner pressure and undergoes material turnover [Holzapfel et al., 2012]. It contains multiple stress-bearing con- stituents, including isotropic elastin-dominated matrix, circumferentially oriented Smooth Muscle Cell (SMC) and collagen-dominated families of fibers with different orientations 7 [Zeinali-Davarani et al., 2011]. In addition, a method of constrained mixture model (CMM) is utilized to homogenize constituents with different constituent-specific deposition pre-stretch. The CMM method allows us to perform finite element simulations on modeling evolution of arterial segments, given as mixtures where each constituent has with distinct residual stress. As a result, there is no need to find a common arterial configuration at zero blood pressure, knowing a fact that the configuration is not usually obtained in human clinical applications. Elastin contributes resilience and elasticity to the aortic tissue, but it degenerates over time and is irreplaceable. The degeneration in elastin causes a localized dilation of the aorta, leading to the weakening of the aortic wall as well as the increase of aortic diameter and wall stress [Ghorpade and Baxter, 1996, He and Roach, 1994, Rizzo et al., 1989]. In addition, the collagen fiber family is suggested to be an important material in supporting the main aortic wall [Choke et al., 2005]. It is continuously removed and produced throughout human life in response to the change of aortic wall stress or strain, thus changing the mechanical properties of the aorta. Given the initial degradation of elastin, the aortic wall stress increases, causing a faster accumulation rate of collagen fiber, which plays a compensation role for the loss of elastin. Therefore, in order to capture the long-term adaptation of blood vessels, we further assume that some of the collagen fiber and SMC continuously generate and degrade in order to restoring toward a value of homeostatic stress [Baek et al., 2006]. Along with the assumption that an irregular elastin degradation initially induces the local dilation of AAA, the proposed G&R model generates long-term changes of the AAA geometrical features, which have been tested for patient-specific simulations [Zeinali-Davarani and Baek, 2012, Zeinali-Davarani et al., 2011]. In this section, we introduce the vascular G&R models by describing its mathematical for- mulations. Next, a pre-step toward the surrogate model is introduced to translate the G&R theories by implementing two models of AAA with different fidelities, where a high-fidelity model is implemented on FEniCS [Alnæs et al., 2015] and a low-fidelity model is imple- 8 mented on MATLAB. A verification between two-fidelity models are provided at the end of this section. 2.2 Mathematical modeling of vascular G&R 2.2.1 Kinematics and configurations Consider a blood vessel as a mixture consisting multiple materials which undergo constant changes including deformation, production and degradation. In CMM approach, we assume that each newly produced constitute is deposited onto the existing vascular tissues with a “deposition pre-stretch” at a homeostatic value. In particular, at the intermediate time τ during the deposition (τ ∈ [0, t]), the natural stress-free configuration of the newly produced k th constituent is denoted by κkn (τ ), and the overall constrained mixture is denoted by κτ . Accordingly, we use a deformation gradient Gk (τ ) to map from the natural configuration κkn (τ ) to the current in-vivo configuration κτ of the k th constituent at time τ . All important configurations and mappings are described in the schematic drawing, Figure 2.1 [Baek et al., 2006, Farsad et al., 2015, Zeinali-Davarani et al., 2011]. For the computational purpose, a reference configuration κR is introduced to provide con- venient mappings among multiple configurations. Note that it is not necessary to prescribe a stress-free reference configuration κR due to the CMM. In this study, the κ0 , the current in vivo configuration at time = 0, is selected to be κR . The mapping from κR to κτ is de- noted by F(τ ). Notice that constituents turnover as time passes, which usually leads to the changes of overall mass at the current configuration, κt . This implies that the mass changes correspond to κR when its evolution is traced over the G&R duration time τ , while κR is fixed in space. 9 Figure 2.1: Important configurations of vascular G&R, including the current in-vivo con- figurations κτ , the computational reference configuration κR , and the natural stress-free configuration of the newly produced constituent κkn (τ ), in which time τ ∈ [0, t]. The de- formation gradient F(τ ) maps from κR to κτ . Gk (τ ) denotes the deposition pre-stretch mapping from κkn (τ ) to κτ . As aforementioned, the CMM theory assumes that multiple constituents compose into a continuous mixture and deform together. More specifically, a single particle of the mixture is contributed by multiple characteristics, such as all types of constituents with different mechanical properties and associated mass. Now we can derive the Fkn(τ ) (t), which represents the deformation gradient for each constituent that maps from its natural configuration at the produced time, i.e., κkn (τ ), to the configuration of mixture at the current time t, i.e., κt . As illustrated in Figure 2.1, Fkn(τ ) (t) is given by Fkn(τ ) (t) = F(t)F−1 (τ )Gk (τ ). (2.1) In the G&R theory, the arterial wall is considered as a homogenized mixture of three stress- bearing constituents, including elastin, collagen fiber and SMC. For convenience, we use 10 subscript i = e, 1, ..., k, m to represent different constituents. In particular, the elastin is denoted by e; the four collagen fiber families are denoted by k = 1, ..., 4; the SMC is denoted by m. Collagen fiber families and SMC are continuously produced and removed over time. Following Figure 2.1, the newly synthesized collagen fibers and SMC are deposited into the tissue with pre-stretch. For simplicity, we assume that all new-produced fibers share the same homeostatic value pre-stretch: Gkh for collagen fiber and Gm h for SMC. For collagen, the direction of k th fiber is represented by a unit vector mk (τ ). The corre- sponding unit vector in the reference configuration κR can be inversely found by F(τ )−1 mk (τ ) Mk (τ ) = . (2.2) |F(τ )−1 mk (τ )| The stretch of collagen fiber (synthesized at time = τ ) caused by the overall deformation of blood vessel is expressed as p λk (t) = F(t)Mk (τ ) · F(t)Mk (τ ). (2.3) Thus, the stretch of collagen fiber maps from the natural stress-free configuration to the current configuration given by λk (t) λkn(τ ) (t) = Gkh , (2.4) λk (τ ) where λk (t) and λk (τ ) denote stretch of collagen fiber maps from the natural to the current configuration at time = t and τ , respectively. The stretch of SMC can be derived in a similar way. Given Equation 2.3, the stretch of SMC that maps from the natural configuration to the current configuration is provided by λm (t) λm n(τ ) (t) = Gh m , (2.5) λm (τ ) in which the SMC is assumed to be distributed in the circumferential direction of the blood 11 vessel. As aforementioned, elastin is another main stress-bearing constituent in the blood vessel. It is produced at a period of prenatal stage and degrades slowly over time compared to the turnover rate of collagen fiber, so it is difficult to identify the deposition pre-stretch of elastin [Langille, 1996]. Thus, a deformation gradient G̃e = F−1 (τ )Geh is defined to represent the mapping between the natural to the reference configuration of elastin. Accordingly, the deformation gradient that maps the elastin from natural to the current configuration is provided as Fen (t) = F(t)G̃e , (2.6) where G̃e is specified by    Ge1 0 0  e   G̃ =  Ge2 0 . (2.7)  0    1 0 0 Ge1 ∗Ge2 2.2.2 Constitutive laws and constituent turnover With the kinematic formulations of constituents, the strain energy for the elastin, collagen fiber families (k =1,...,4) and passive SMC per unit mass are given as   e c1 1 Ψ (Cen (t)) = e e C[11] + C[22] + e e e 2 −3 , (2.8) 2 C[11] C[22] − C[12] k k c2 n k 2 2 o Ψ (λn(τ ) (t)) = exp[c3 (λn(τ ) (t) − 1) ] − 1 , (2.9) 4c3 c4 n m 2 2 o Ψm (λm n(τ ) (t)) = exp[c (λ 5 n(τ ) (t) − 1) ] − 1 , (2.10) 4c5 e e e where C[11] , C[22] and C[12] are components of Cen (t) = (Fen (t))T Fen (t) which represents the Green tensor of elastin mapping from its stress-free natural configuration to the current 12 configuration. Additionally, due to the constant turnover of vascular tissue, the mass per area with respect to the reference configuration, i.e., M (t), can be given by X X Z t  M (t) = i M (t) = MRi (0)Qi (t) + miR (τ )q i (t, τ )dτ (2.11) i i 0 where Qi (t) and q i (t, τ ) represent the survival functions, which are the survival fraction (the part of material which has not degrade) of the initial mass and the mass produced at time τ , respectively. MRi (0) represents the initial mass of ith constituent at reference configuration, and miR (τ ) represents the stress-mediated mass production rate, which will be further explained in Section 2.2.3. Now, given the mass per area, we can give the share of thickness for each constituent with respect to both the reference configuration, H i (t), and the current configuration, hi (t), by M i (t) M i (t) H i (t) = , hi (t) = , (2.12) ρ Jρ where J = detF(t) and ρ denotes the volume density of the vascular tissue. The total thickness is then given by X X H(t) = H i (t), h(t) = hi (t). (2.13) i i Furthermore, the survival function q i (t, τ ) is specified by    exp(− τ k i (τ̃ )dτ̃ ) t − τ ≤ ai  R   t q max q i (t, τ ) = , (2.14)   0 t − τ > aimax   where kqi (τ̃ ) and aimax denote the rate of removal and maximum lifespan to ith constituent, respectively. 13 For each constituent, the energy per unit reference area of AAA is expressed as the summa- tion of all strain energies contributed by various survival masses, i.e., we (t) = M e (t)Qe (t)Ψe (Cen (t)), (2.15) Z t k k k k w (t) = M (0)Q (t)Ψ (Ckn(0) (t)) + mk (τ )q k (t, τ )Ψk (Ckn(τ ) (t))dτ, (2.16) 0 Z t m m m m w (t) = M (0)Q (t)Ψ (Cmn(0) (t)) + mm (τ )q m (t, τ )Ψm (Cm n(τ ) (t))dτ, (2.17) 0 where Ψ(Cin(τ ) (t)), i = e, k, m, represents the stored energy of constituent i synthesized at time τ , and the Green tensor Cin(τ ) (t) = {Fkn(τ ) (t)}T Fkn(τ ) (t). Therefore, the total strain energy per unit reference area can be given by X w(t) = we (t) + wk (t) + wm (t). (2.18) k 2.2.3 Homeostasis and stress-mediated adaption The Cauchy stress of the mixture is provided by X t= ti , (2.19) i where for each constituent i, we have 2 ∂wi ti = F(t) F(t)T . (2.20) J ∂C(t) Then, the average value of the stress for k th collagen fiber that is measured by a scalar value is provided by ||tc (t)mk (t)|| σ k (t) = , (2.21) hc (t) 14 where tc (t) and hc (t) are the summation of Cauchy stress tensor and the thickness for all collagen fiber families combined. The scalar stress of SMC can be calculated in a similar way. In order to capture the stress-mediated adaption, we further assume that the mass production rate of collagen fiber and SMC are proportional to the difference between scalar stress of fiber and a homeostatic stress σhi [Baek et al., 2006], i.e., MRc (t) k k mkR (t) = (K (σ (t) − σhc ) + mkbasal ), (2.22) MRc (0) g MRm (t) m m mmR (t) = (K (σ (t) − σhm ) + mm basal ), (2.23) MRm (0) g where MRc (t) is the total mass density of the collagen fiber family in reference configuration, and Kgi is a scalar parameter which controls the magnitude of the stress-mediated mass production rate. Accordingly, a larger Kgi implies that the blood vessel is able to produce more collagen fiber to maintain the stability of mechanical properties under the elevated vascular stress from elastin degeneration. Therefore, Kgi plays a decisive role in controlling the self-repairing and evolutionary process of blood vessel. 2.2.4 Elastin degradation and smooth muscle active tone The degeneration in elastin causes a localized dilation of the aorta, leading to the weakening of the aortic wall as well as the increase of aortic diameter and wall stress. In this study, the degradation in elastin is prescribed by a spatial damage function which is axisymmetric to the vascular centerline. Specifically, the ratio of elastin degradation at different spatial locations on the vascular centerline (s) is specified by a Gaussian form function: " # 2 (s − µdmg ) dinit = exp − 2 , (2.24) 2σdmg 15 where µdmg represents the center of the Gaussian function with the maximum elastin degra- dation; σdmg defines the standard deviation of the Gaussian function, which controls the area of degraded elastin. µdmg and σdmg affect the initial loss of elastin, which could further lead to the changes of stress-stretch and geometrical state of the AAA, causing various geometrical features of the AAA. The initial damage function of elastin in Equation 2.24 initializes the development of AAA in its early stage. In addition, Wilson et al. [2012], Zeinali-Davarani et al. [2011] introduced time-dependent degradation elastin, which assumes additional removal to elastin caused by the damage of AAA when vascular tissue is under large stretch. This ‘stretch-induced’ damage of elastin is employed to capture the increase of expansion rate in the later stage of AAA development. The over-stretch damage function of elastin is provided by       1 (I1e − 3) ≥ 4.48         π(7.48−I1e ) ddmg = 1 − sin 2.96 3 ≤ 1(I1e − 3) < 4.48 , (2.25)       0 (I1e − 3) < 3     where I1e is the first invariant of the Green tensor Cen . The overall survival elastin is expressed by (1 − dinit )(1 − ddmg ). Accordingly, the mass of elastin M e during G&R process in the reference configuration is provided by M e (t) = M e (0)(1 − dinit )(1 − ddmg ), (2.26) where M e (0) represents the initial mass of elastin. Furthermore, due to the intrinsic in- terrelationship between elastin and SMC [Karnik et al., 2003, Li et al., 1998], the ratio of removed SMC changes proportional to the ratio of degraded elastin. We take this effect into account by reducing the associated energy and mass of SMC proportionally with the elastin degradation. 16 It is not enough to prescribe the mechanical properties of SMC by using its passive strain energy. The active tone of SMC also plays an important role in adjusting the flow and pressure under different environmental stimulus. Thus, an additional function is provided by Baek et al. [2007], Zulliger [2004] to represent the stain energy of SMC active tone, i.e., 1 (λM − λ2 (t))3   S Ψmact (t) = λ2 (t) + , (2.27) ρ 3 (λM − λ0 )2 where λM and λ0 represent the maximum and minimum value of SMC active tone, and S denotes the maximum stress caused by active tone. Thus, considering the SMC active tone, the total strain energy of the mixture per unit area with respect to reference configuration is given by X w(t) = we (t) + wk (t) + wm (t) + M m (t)Ψm act (t). (2.28) k 2.3 High-fidelity 3D membrane G&R model The proposed vascular G&R has been demonstrated by computational models such as a 2D asymmetrical model [Baek et al., 2006] and a 3D membrane model [Zeinali-Davarani and Baek, 2012, Zeinali-Davarani et al., 2011] on MATLAB. However, they are either too simple or too computational costly for practical uses. In this section, a 3D membrane G&R model is established using FEniCS [Alnæs et al., 2015] (a finite element method software implemented on Python) to enable fast G&R simulations. 17 2.3.1 Governing equations In this section, we use the principle of virtual work [Kyriacou et al., 1996] to formulate the weak form of the proposed strained energy, i.e., Z Z δI = δwdA − P n · δxda = 0, (2.29) S s where P n denotes the inner pressure vector applied to a vascular wall; S and s correspond to the surface of a blood vessel in reference and current configurations, respectively; δx denotes the virtual changes in position. Note that the energy form changes with the current time t, so Equation (2.29) needs to be solved at every time step. Next, the current location is replaced by finite element approximation, i.e., δx = Φδxp , (2.30) where xp is the nodal vector for the current position and Φ is the shape function matrix. The governing equation defined on nodal points of local element is then formulated by Z   ∂w ∂Cαβ {F}eP = − P̃i ΦiP dA = 0, (2.31) S ∂Cαβ ∂xpP where xpP is the index notation of current position vector xp , and P̃i is given by P̃ = P JF−T N, (2.32) where N is the normal vector of element in the reference configuration. To solve Equation (2.31), we use the Newton-Rahpson method which requires a tangent 18 matrix defined by ∂ 2w ∂Cαβ ∂Cγω ∂w ∂ 2 Cαβ h ie h ∂F ie Z   K = = − P̃i,Q ΦiP dA = 0. (2.33) PQ ∂xp PQ Se ∂Cαβ ∂Cγω ∂xpP ∂xpQ ∂Cαβ ∂xpP ∂xpQ Note that here we use a complex notation system in which i, j, k, l, m are three-dimensional indexes of the global coordinate system; α, β, γ, ω are two-dimensional indexes of the local coordinate system; A, B, P, Q are indexes of np -dimensional array in a local element (np is the degree of freedom in each element). 2.3.2 Membrane element and local coordinate system Given the global coordinate X = {X1 , X2 , X3 } defined on the base vector {E1 , E2 , E3 } in the reference configuration, we can access to the nodal coordinates of each linear triangular membrane element, i.e., {X(1) , X(2) , X(3) }, from a meshed geometry. In addition, the global coordinates on the deformed geometry are defined on the same Cartesian coordinate system, and represented by x = {x1 , x2 , x3 } in the current configuration. The associated nodal points are given by {x(1) , x(2) , x(3) }. In addition to the general finite element formulation, the fiber orientations of anisotropic constituents should be specified in each local element. In particular, the centroid Xc and the outward normal vector Ne can be directly computed from the nodal coordinates of elements; the local base Ee1 is prescribed by the axial direction of the blood vessel; local base Ee2 , which represents the circumferential direction of the fiber, is computed by Ee2 = Ne × Ee1 . (2.34) Given the basis of the Cartesian coordinate of each local element {Ee1 , Ee2 , Ne }, the local coordinate in the reference and current coordinates are denoted by Xe = {X1e , X2e , X3e } and 19 Figure 2.2: X = {X1 , X2 , X3 } are global Cartesian coordinates defined on base vector {E1 , E2 , E3 } in reference configuration. {Ee1 , Ee2 , Ne } are transformed base vectors that define the local coordinates. The centroid Xc and nodal coordinates of elements {X(1) , X(2) , X(3) } are defined by global coordinates. The associated nodal coordinates in current configuration of G&R process are provided by {x(1) , x(2) , x(3) }. xe = {xe1 , xe2 , xe3 }, respectively. 2.3.3 Finite element formulations In FEniCS, the shape functions {φi } of a linear triangular element is automatically prescribed such that the values of shape functions on the local element follow (j) (j) (j) φi (X1 , X2 , X3 ) = δij . (2.35) Accordingly, the current location of any point of the geometry is approximated by x = Φxp . (2.36) 20 where xp denotes the nine-dimensional global nodal coordinates of the linear triangular element (1) (1) (1) (2) (2) (2) (3) (3) (3) xp = {x1 , x2 , x3 , x1 , x2 , x3 , x1 , x2 , x3 }T , (2.37) and Φ denotes the 9 × 3 shape function matrix   1 2 3  φ 0 0 φ 0 0 φ 0 0    Φ(X1 , X2 , X3 ) =   0 φ1 0 0 φ2 0 0 φ3 0 . (2.38)   0 0 φ1 0 0 φ2 0 0 φ3 Now, given the global coordinates in the current configuration, i.e., x, enabled by FEniCS, we need to formulate the local coordinates and all forms in Equation (2.31) and Equation (2.33) to find the solution. Given the base vectors of both global and local coordinate systems, i.e., {E1 , E2 , E3 } and {Ee1 , Ee2 , Ee3 } (let Ee3 = Ne for convenience ), the transformation formula is defined as Xe = Q(X − Xc ), (2.39) xe = Q(x − Xc ), (2.40) where the transformation matrix Q is defined by    Ee1 · E1 Ee1 · E2 Ee1 · E3    Q= e e  E2 · E1 E2 · E2 E2 · E3 e .  (2.41)   Ee3 · E1 Ee3 · E2 Ee3 · E3 Accordingly, the local coordinate can be presented by xe = Q(x − Xc ) = Q(Φxp − Xc ), (2.42) 21 i.e., xei = Qij (ΦjA xpA − Xjc ). (2.43) The spatial derivative of location position with respect to the local basis becomes ∂ xei,α = [Qij (ΦjA xpA − XjC )] ∂Xαe ∂ ∂Xk = [Qij (ΦjA xpA − XjC )] ∂Xk ∂Xαe (2.44) ∂Xk = Qij ΦjA,k xpA ∂Xαe = Qij ΦjA,k xpA Qαk . Next, the Gauchy green tensor is given by ∂xe ∂xe   C= · Eeα ⊗ Eeβ , (2.45) ∂sα ∂sβ which can be computed using FE formulation, and shown by Cαβ = xei,α xei,β = ΦeiA,α xpA ΦeiB,β xpB = Qij ΦjA,k xpA Qαk Qil ΦlB,m xpB Qβm (2.46) = δjl ΦjA,k xpA Qαk ΦlB,m xpB Qβm = ΦjA,k xpA Qαk ΦjB,m xpB Qβm . The Gauchy-Green tensor can also be simplified into a matrix notation that is directly implemented in FEniCS, i.e., C = Q̄(∇x)T (∇x)Q̄T . (2.47) Note that the Gauchy-Green tensor C here is a 2 × 2 matrix defined on local coordinate, 22 and Q̄ is a 2 × 3 matrix defined by   e e e  E1 · E1 E1 · E2 E1 · E3  Q̄ =  . (2.48) e e e E2 · E1 E2 · E2 E2 · E3 ∂Cαβ ∂ 2 Cαβ Next, we derive ∂xpP and ∂xpP ∂xpQ defined in Equation (2.31) and Equation (2.33) to enable the the solver. ∂Cαβ ∂ p p p = p (ΦjA,k xA Qαk ΦjB,m xB Qβm ) ∂xP ∂xP = δAP ΦjA,k Qαk ΦjB,m xpB Qβm + δBP ΦjA,k xpA Qαk ΦjB,m Qβm (2.49) = ΦjP,k Qαk ΦjB,m xpB Qβm + ΦjA,k xpA Qαk ΦjP,m Qβm = Qαk (ΦjP,k ΦjB,m xpB + ΦjA,k xpA ΦjP,m )Qβm , ∂ 2 Cαβ ∂ p p = (Qαk (ΦjP,k ΦjB,m xpB + ΦjA,k xpA ΦjP,m )Qβm ) ∂xP ∂xQ ∂xpQ ∂ (2.50) = (Qαk (δBQ ΦjP,k ΦjB,m + δAQ ΦjA,k ΦjP,m )Qβm ) ∂xpQ = Qαk (ΦjP,k ΦjQ,m + ΦjQ,k ΦjP,m )Qβm . The matrix notation leads to a concise form by ∂C  T T  = Q̄ (∇Φ (P ) ) (∇x) + (∇x)) (∇Φ (P ) ) Q̄T , (2.51) ∂xpP ∂ 2C  T T  = Q̄ (∇Φ (P ) ) (∇Φ (Q) )) + (∇Φ (Q) )) (∇Φ (P ) ) Q̄T , (2.52) ∂xpP ∂xpQ where Φ(P ) is the P th column of 3 × 9 shape matrix Φ. 23 Figure 2.3: A verification between two 3D membrane G&R simulations implemented on FEniCS and MATLAB, respectively. The same simulation setup and input parameters are used in two simulations. 2.3.4 Code verification To verify the FEniCS code of 3D membrane AAA G&R model, another MATLAB-based 3D membrane G&R simulation is introduced. A comparison between two simulations with the same setup and input parameters are provided in Figure 2.3. The verification result shows that two AAAs share the same size of geometries and similar stress distributions. 2.4 Low-fidelity axisymmetric G&R model and verifi- cation The axisymmetric G&R model uses a cylindrical polar coordinate system, wherein the refer- ence and current coordinates are defined by X = (Xz , Xr ) and x = (xz , xr ), respectively. In addition, the initial shape of artery is given by Xr = R(Xz ) in the reference configuration. 24 x= (xz, xr) E1 x = (x1, x2, x2) Local element E2 Local element z r 3D G&R model Axisymmetric G&R model Figure 2.4: Two-fidelity G&R FEM models. The 3D G&R model is built on a three dimen- sional geometry of artery with linear membrane elements. The axisymmetric G&R model is built on an idealized axisymmetric geometry with quadratic membrane elements. The comparison between two G&R models is illustrated in Fig. 2.4. A verification results of AAA growth simulated by the high-fidelity and low-fidelity vascular G&R models are shown in Fig. 2.5. The same input parameters are used in both models. Accordingly, for axisymmetric G&R model, the governing equation is      2π L ( ∂w Φi + ∂w 0 p P x0z xr Φi dXz R Φ )R 1 + (R0 )2 −  0       0 ∂xr ∂x0r i F= = , (2.53)  2π L ( ∂w Φ + ∂w 0 p Φ )R 1 + (R0 )2 + P x0r xr Φi dXz  R  i    0   0 ∂xz ∂x0 i z where L is the max value of Xz . 25 (a) Maximum Diameter against time (t) (b) Diameter against length(s) Figure 2.5: Verification results of AAA growth between the FEniCS-based high-fidelity model (represented by ‘FEniCS’) and the MATLAB-based low-fidelity model (represented by ‘Mat- lab’). The maximum diameter (max. d) vs. time (t) for both simulations are plotted in Figure a. The diameter curves (d) along the length of arterial centerline (s) for both simu- lations at t = 2000 days and t = 3000 days are plotted in Figure b. Note that the crosses and circles in Figure b represent the AAA diameters at nodal locations. 26 Chapter 3 PHYSICS-BASED SURROGATE MODELS 3.1 Introduction Patient-specific applications of some computational vascular models, such as vascular G&R models, are too computationally expensive for practice in clinical treatment and decision making. This is mainly because the patient-specific modeling needs a series of time-consuming operations, such as image segmentation, meshing, patient-specific parameter estimation and model validation. In particular, the patient-specific parameter estimation1 is particularly computationally expensive because it is considered as a numerical inverse problem that re- quires extensive forward simulations to solve [Marsden, 2013]. For an instance, it sometimes takes days and even months to estimate the patient-specific parameters for a single patient. In contrast, the clinical practice requires the personalized suggestion to be provided within the clinical relevant time frame (hours) [Taylor et al., 1999]. Thus, despite the existence of many significant patient-specific vascular models [Caroli et al., 2013, Colombo et al., 2019, 27 Kuhl et al., 2007, Long et al., 1998, Zambrano et al., 2017, Zeinali-Davarani and Baek, 2012], direct access to rapid clinical recommendations remains challenging. Therefore, there is a pressing need to design an efficient patient-specific simulation and prediction framework to aid clinical treatments of vascular diseases. To achieve this goal, some studies employed reduced-order models for parameter searching [Itu et al., 2015, Spilker and Taylor, 2010, Zhang et al., 2019], establishing efficient opti- mization processes, but with reduced accuracy. In particular, the Probabilistic Collocation Method (PCM), a computationally efficient approximation approach, interpolates compu- tationally expensive simulations sampled on collocation points and predicts intermediate simulation results. PCM has been broadly used in climate research, energy system, hemo- dynamics. In Section 3.2, we utilize PCM as a fast surrogate model to augment the dataset simulated by AAA G&R model. In addition, the Kriging-based surrogate model, a cheap-to-evaluate machine learning model that approximates the original model, has been broadly applied for efficient optimization [Dubourg et al., 2011, Li and Yang, 2019, Matheron, 1963, Sacks et al., 1989]. In the field of arterial biomechanics, the Kriging method has been successfully demonstrated by providing surrogates of physical models of cardiovascular systems in various studies including the identification of design parameters of bypass graft [Sankaran and Marsden, 2010], scaffold [Szafron et al., 2019] and baffle [Yang et al., 2010], homeostatic parameters in arterial growth and remodeling [Sankaran et al., 2013], and geometrical parameters of arterial branches [Marsden et al., 2008]. Although advanced approximation techniques such as Kriging accelerate the time-consuming cardiovascular simulations, the creating of surrogate models needs a large number of training data, i.e., a large set of forward vascular simulations. The generation of training data may 1 In this study, the ‘patient-specific parameter estimation’ particularly refers to the process of estimat- ing some parameters in a vascular model by minimizing the discrepancy between simulation outputs and personalized clinical data. 28 still be overwhelmed expensive due to the complexity of the vascular model. Fortunately, cokriging, i.e., a multi-variable extension of kriging, can efficiently reduce the amount of expensive sample data required for model training. In particular, cokriging is the multivariate Gaussian Process (GP) regression that treats training samples as deterministic measurements in the parameter space, and predicts the intermediate value and related error at any arbitrary location [Fernandez et al., 2017, Kennedy and O’Hagan, 1998, Sacks et al., 1989]. Cokriging provides an efficient metamodel approach of combining the simulations of High-Fidelity Model (HFM) with Low-Fidelity Model (LFM) [Forrester et al., 2007]. In Section 3.3, we consider HFM as the highly expensive computational model with an adequate accuracy, which can be remedied by a large amount of efficient simulations from the less accurate LFM via cokriging. By introducing HFM, LFM, and cokriging, we can construct an efficient Multi-fidelity Surrogate (MFS) that plays as an important middle stage connecting the vascular G&R simulations to the personalized parameter estimation and prediction. In this chapter, the theories of PCM, kriging and cokriging are described to enable useful surrogate models, which further lead to more complex operations such as patient-specific parameter estimation (Chapter 4) and prediction (Chapter 5 and Chapter 6). 3.2 Probabilistic collocation method 3.2.1 Deterministic input Assume the physics model takes a group of input parameters γ = {γ [1] , γ [2] , γ [3] , ...} and generates simulation output y = η(γ), (3.1) 29 where η(.) represents the physics-based computational model. To reduce the computational cost, we shall approximate η(.) by utilizing a set of N basis functions {gi (γ)}, with i = 1, · · · , N , such that XN ŷ = βi gi (γ), (3.2) i=0 where βi represents the regression coefficient. Given a set of functions {gi (γ)}, the regression coefficients {βi } can be solved as follow. The residual between the truth and the approximation is defined as R({βi }, γ) = ŷ(γ) − y(γ). (3.3) By applying the ordinary least squares estimation to (3.2), the optimal set of coefficients βi is formulated in Z hgi (γ), R({βi }, γ)i = gi (γ)R({βi }, γ)dγ = 0, (3.4) γ where i = 1, · · · , N and h., .i represents the dot product between two deterministic functions. (3.4) can be solved by the idea of Gaussian quadrature [Stroud and Secrest, 1966], which approximates the integral as Z XN R({βi }, γ)gi (γ)dγ ' vj R({βi }, γ̃j )gi (γ̃j ) = 0, (3.5) γ j=1 where vj are weights and γ̃j are abscissas, respectively. If the weights and the basis functions Q are chosen such that i,j vj gi (γ̃j ) > 0 for all i and j, the summation in (3.5) can be further approximated as R({βi }, γ̃j ) = 0, j = 0, · · · , N. (3.6) Note that the quadrature points γ̃j are also the collocation points. (3.6) can be used to find the coefficients {βi } by running tfhe model at N + 1 different collocation points and solving 30 a system of N + 1 equations. 3.2.2 Stochastic input Suppose that the input γ is a random vector with a known probability density function (PDF) π(γ), (3.4) can be transformed into the probability space as Z π(γ)R({βi }, γ)gi (γ)dγ = 0. (3.7) γ Similarly, with the proper choice of vj and gi (γ), (3.6) becomes π(γ)R({βi }, γ̃j ) = 0, (3.8) where j = 0, · · · , N . Since the PDF function π(γ) is always positive, (3.6) can be used to find the coefficients in the stochastic case. 3.2.3 Selection of base functions and collocation points Theorem 3.1 Consider a quadrature formula: Z X N W (γ)F (γ)dγ ' wj F (γj ), (3.9) γ j=1 where wj are weights and γj are abscissas. Given a weight function W (γ) = γ α (1 − γ)β , an optimal choice of N quadrature points can be found based on the correct integration using the highest order of the polynomial expansion of F (γ). The optimal quadrature points are the zeros of the polynomial of degree (N + 1), i.e., PN +1 (γ), that satisfy the orthogonality 31 condition Z W (γ)γ j PN +1 (γ)dγ = 0, for j = 0, · · · , N. (3.10) γ The detailed proof of Theorem 3.2.3 is provided in Chap 3 of [Villadsen and Michelsen, 1978]. In short, the choice of collocation points as the roots of the next order orthogonal polynomial lets the collocation method approximation close to Galerkin’s method. As a result, it outperforms other methods of weighted residual (MWR) [Villadsen and Michelsen, 1978]. Corollary 3.2 Consider the same quadrature formula in Theorem 3.2.3, and a set of N + 1 orthogonal polynomial functions {gi (γ)}, if the set of functions satisfy the condition Z π(γ)gi (γ)gN +1 (γ)dγ = 0, i = 1, · · · , N, (3.11) γ they also satisfy condition (3.10), and the zeros of gN +1 (γ) are the optimal quadrature points. Proof of Corollary 3.2: We illustrate the proof for N = 2, then the proof for an arbitrary N is straightforward. Let g0 = 1, g1 = ax + b, g2 = cx2 + dx + e. 32 So, applying 3.11, we have Z Z π(x)g0 (x)g3 (x)dx = π(x)g3 (x)dx = 0, x x Z Z π(x)g1 (x)g3 (x)dx = π(x)(ax + b)g3 (x)dx = 0, x x Z Z π(x)g2 (x)g3 (x)dx = π(x)(cx2 + dx + e)g3 (x)dx = 0. x x Rearrange the terms, we have an equivalent system of equations Z Z (1 + b + e) π(x)g3 (x)dx = 0 → π(x)g3 (x)dx = 0, x Z Z x (a + d) π(x)xg3 (x)dx = 0 → π(x)xg3 (x)dx = 0, Z x Z x (c) π(x)x2 gx (x)dx = 0 → π(x)x2 g3 (x)dx = 0. x x which is the condition 3.10. If we choose the weight function to be the PDF of γ, i.e., W (γ) = π(γ), (3.11) can be used to generate the set {gi (γ)} in a recursive manner as follows. In practice, we define the initial conditions g−1 = 0, g0 = 1, and the orthogonal polynomials can be obtained recursively by solving the equations Z π(γ)gi (γ)gi+1 (γ)dγ = 0, i = 1, · · · , N. (3.12) γ However, solving for high order polynomials (3.11) is time-consuming and error-prone. Thus, we use Favard theorem to compute the set of basis functions more efficiently. 33 Theorem 3.2 (Favard Theorem) If a sequence of polynomials {Pi (γ)}, where i = 1, · · · , N , satisfies the recurrence relation Pi (γ) = (γ − ci )Pi−1 (γ) − di Pi−2 (γ), (3.13) P−1 = 0, P0 = 1, where ci and di are real numbers, then {Pi (γ)} are orthogonal polynomials. Algorithm 1 The Probabilistic Collocation Method (PCM) model Part 1: Compute collocation points 1: Initialize g−1 (γ) = G−1 (γ) = 0 and g0 (γ) = G0 (γ) = 1. 2: for i = 1, · · · , N do p 3: Gi (γ) = γgi−1 (γ) − hγgi−1 (γ), gi−1 (γ)igi−1 (γ) − hGi−1 (γ), Gi−1 (γ)igi−2 (γ) Gi (γ) 4: gi (γ) = √ hGi (γ),Gi (γ)i 5: end for 6: Find the zeros of gN +1 (γ) = 0 as N collocation points. 7: Repeat steps 1-5 for other random variables. We denote three random variables2 and three sets of basis functions as [1] [2] [3] [1] [2] [3] {γi , γi , γi } and {gj (γ), gj (γ), gj (γ)}, respectively, where i = 1, · · · , N and j = 0, · · · , N − 1 8: Create a permutation of three groups of collocation points of 3 random variables, i.e., (γi[1] , γj[2] , γk[3] ) where i = 1, · · · , N , j = 1, · · · , N , and k = 1, · · · , N . Part 2: Run the computational code at the collocation points 1: for i = 1, · · · , N , j = 1, · · · , N , k = 1, · · · , N do [1] [2] [3] 2: Run η(γi , γj , γk ). 3: end for Part 3: Compute the coefficients β  [1] [2] [3]  η(γ1 , γ1 , γ1 ) 1: Concatenate the computational outcomes: y =   ..  .  .  [1] [2] [3] η(γN , γN , γN )  [1] [1] [2] [2] [3] [3] [1] [1] [2] [2] [3] [3]  gN −1 (γ1 )gN −1 (γ1 )gN −1 (γ1 ) · · · g0 (γ1 )g0 (γ1 )g0 (γ1 ) 2: Arrange the matrix K =   .. .. ..  .  . . .  [1] [1] [2] [2] [3] [3] [1] [1] [2] [2] [3] [3] gN −1 (γN )gN −1 (γN )gN −1 (γN ) · · · g0 (γN )g0 (γN )g0 (γN ) 3: Compute the coefficients: β = K −1 y. (3.14) Part 4: Approximate the computational code 1: For any new set of random variable (γ∗[1] , γ∗[2] , γ∗[3] ), the outcome can be approximated by: N −1 N −1 N −1 [1] [2] [3] [1] [1] [2] [2] [3] [3] X X X η ∗ (γ∗ , γ∗ , γ∗ ) = β i,j,k gi (γ∗ )gj (γ∗ )gk (γ∗ ). (3.15) i=0 j=0 k=0 p In this study, we utilized ci = hγgi−1 , gi−1 i and di = hgi−1 , gi−1 i by referring the work 2 For the sake of notational simplicity, we denote {Kg , σd , µd } as {γ [1] , γ [2] , γ [3] } in Algorithm 1. 34 of [Zhou et al., 2014]. The overall PCM algorithm is shown in Algorithm 1. Furthermore, though the PCM has been widely used to approximate a univariate prediction target, i.e., a scalar value y, it is extended to be a multivariate approximation in this study. 3.3 Multi-fidelity surrogate In general, high-fidelity computational vascular simulations are computationally expensive due to high complexity of the arterial adaptation mechanism. For examples, using a 3.3 GHz CPU and a 64 GB RAM, a 3D G&R simulation of AAA implemented on MATLAB takes around two days [Zeinali-Davarani and Baek, 2012]; a 3D hemodynamic simulation of a blood vessel in a cardiovascular cycle takes approximately 6 hours on a 16 cores 2.67 GHs GPU. As those examples shown, even the computational time of a single forward vascular simulation is significantly larger than the clinical relevant time (1-2 hours). Consequently, it becomes computationally impractical to perform more time-consuming applications, such as personalized parameter estimation which could be considered as an inverse problem that requires extensive forward vascular simulations to solve. To address this issue, we use a surrogate model to approximate computationally expensive physics-based simulations, which significantly reduces computational time. A good option of creating surrogate model is to use kriging, i.e., Gaussian Process (GP) regression, which treats a set of sample data as deterministic measurements, and predicts the intermediate value and the associated variance at an arbitrary location in the parameter space [Sacks et al., 1989]. The kriging method, however, is limited as a surrogate model because it can only approximate a single computational model, whereas we can access to multiple vascular models with different fidelities. Consequently, we employ cokriging, i.e., a multivariate extension of the single variate kriging3 , to combine the information of multi- 3 The kriging is particularly referred to the univariate GP regression in the context of this study. 35 fidelity models, and thereby further improving the computational efficiency. Here, we denote the approximated model established by cokringing as MFS. A cokriging-based MFS can be trained by a sufficient number of computationally cheap but less accurate simulations and a limited number of expensive but accurate simulations, and thus reaching a balance between accuracy and efficiency. We consider the creating of corking- based MFS as a necessary intermediate process between the vascular simulations (Chapter 2) and the patient-specific parameter estimation (Chapter 4). For detailed derivations of cokriging-based MFS, readers can refer to Forrester et al. [2007], Jones et al. [1998], Kennedy and O’Hagan [1998]. 3.3.1 Kriging Kriging, the univariate version of GP regression, is introduced first to illustrate the basic concepts in cokriging. Essentially, kriging is an approximation method which interpolates a set of sample data. Consider a physical process, e.g., a vascular simulation, that generates a target value y = f (x) under an input column vector of a sample x ∈ RNx . In a set of Ne samples, the sample locations and corresponding target values are provided by X = {x(1) , ..., x(Ne ) } and y = {y (1) , ..., y (Ne ) }T , respectively. Kriging treats each sample value y as a realization of a stochastic process Y (x), which provides the best linear unbiased prediction in a regression framework [Sacks et al., 1989]. For ordinary kriging, the stochastic process is prescribed by Y (x) = µ + Z(x), (3.16) which hypothesizes that the Y (x) is modeled as the summation of an unknown constant hyperparameter µ and a GP term Z(x) with zero mean. Z(x) can be interpreted as the local feature of the approximated stochastic process, i.e., the local deviation from the average value µ. Moreover, if the approximated physical process f (x) is continuous and smooth, it 36 is natural to suppose that the covariance of Z(x) is related to distance [Jones et al., 1998]. This idea of covariance is originally proposed in geostatistics to optimize the mining locations [Matheron, 1963]. Accordingly, the covariance of Z(x) is specified by a Gaussian correlation function defined purely on a weighted spatial distance,  X Nx  ∗ 2 ∗ 2 ∗ 2 cov{Z(x), Z(x )} = σ ψ(x, x ) = σ exp − pi (xi − xi ) , (3.17) i=1 where hyperparameters include variance σ and the magnitude of variation in different dimen- sions p = {p1 , ..., pNx }. This covariance function guarantees that the correlation converges to one when x and x∗ come close enough. The correlation function, denoted by ψ in (3.17), is also used for constructing the Nx × Nx correlation matrix   (1) (1) (1) (Nx )  ψ(x , x ) ... ψ(x , x )   .. .. ..  Ψ=  . . . ,  (3.18)   ψ(x(Nx ) , x(1) ) ... ψ(x(Nx ) , x(Nx ) ) and correlation vector function ψ(x) = {ψ(x, x(1) ), ..., ψ(x, x(Nx ) )}T , given a set of sample points X. In short, the stochastic process has been clearly defined, and parameterized by hyperparam- eters, including µ, σ and p. In practice, those hyperparameters are required to be carefully calibrated such that the conditional posterior probability of Y (x) is maximized given the sample data, so the stochastic process can optimally approximate the physics-based model at the sample points. In particular, the estimations of µ and σ are provided by µ = 1T Ψ−1 y(1T Ψ−1 1)−1 , (3.19) σ 2 = (y − µ1T )Ψ−1 (y − µ1)/Nx , (3.20) 37 and the p can be computed by maximizing the concentrate ln-likelihood n 1 − ln σ 2 − ln | det(Ψ)|. (3.21) 2 2 Now, with the hyperparameters calibrated by sample data, the best linear unbiased predic- tion ŷ(x) and its variance σy2 (x) are provided as ŷ(x) = µ + ψ T (x)Ψ−1 (y − µ1), (3.22) (1 − 1T Ψ−1 ψ(x))2   2 2 T −1 σy (x) = σ 1 − ψ (x)Ψ ψ(x) + , (3.23) 1T Ψ−1 1 respectively, where 1 denotes a column vector with Nx ones. Note that the variance is given in the form of standard deviation. 3.3.2 Cokriging Cokriging is also called multivariate GP regression. It has been applied to many problems for which multilevel analysis is available. Consider two physics-based models that simulate for the same process. Models that are computationally expensive but more accurate are labeled as High-fidelity models (HFMs); less expensive but more accurate are labeled as Low-fidelity Models (LFMs). In practice, typical LFMs are cheaper because they are usually simplified from HFMs. For example, compared to HFM, LFM may hold idealized geometry, smaller mesh size, simpler physics assumption, reduced dimensionality, etc. In our study, given an input column vector x ∈ RNx , the HFM produces a computational expensive output ye = fe (x), while the LFM produces a cheap output yc = fc (x). Similar to ordinary Kriging (3.16), cokriging also assumes that, for both HFM and LFM, the approx- imated stochastic process is the sum of a constant and a GP with zero mean. In addition, 38 the cokriging further prescribes the relation of GPs by Ze (x) = ρZc (x) + Zd (x), (3.24) where Ze (x) and Zc (x) denote the GPs of HFM and LFM, respectively; Zd (x) is defined as a GP that is independent of Zc (x) and represents the difference; ρ is a hyperparameter. Moreover, the covariance, correlation function, correlation matrix and correlation vector function of Zc (x) and Zd (x), are defined in the same manner as their counterparts shown in Section 3.3.1. For convenience, variables and parameters in different stochastic processes are represented by different subscripts: e for the HFM related, c for the LFM related, and d for Zd (x) related. Suppose that we have a set of Ne High-fidelity (HF) sample simulations with inputs Xe = (1) (Ne ) (1) (Ne ) T {xe , ..., xe } and outputs ye = {ye , ..., ye } , as well as a set of Nc Low-fidelity (LF) (1) (Nc ) (1) (Nc ) T sample simulations with inputs Xc = {xc , ..., xc } and outputs yc = {yc , ..., yc } . Based on these sample simulations, the hyperparameters {µc , σc2 , pc , µd , σd2 , pd , ρ} in cokriging can be estimated using MLE in a similar way to (3.19), (3.20) and (3.21). The readers could follow Kennedy and O’Hagan [1998] and Forrester et al. [2007] for detailed descriptions. Therefore, given the complete form of covariance matrix   σc2 Ψc (Xc , Xc ) ρσc2 Ψc (Xc , Xe ) Ψ= , (3.25)   ρσc2 Ψc (Xe , Xc ) ρ2 σc2 Ψc (Xe , Xe ) + σd2 Ψd (Xe , Xe ) and the complete covariance vector   ρσc2 ψc (Xc , x) ψ(x) =  , (3.26)   ρ2 σc2 ψc (Xe , x) + σd2 ψd (Xe , x) 39 the cokriging-based best linear unbiased prediction and the associated variance are ŷ(x) = µ + ψ T (x)Ψ−1 (y − µ1), (3.27) (1 − 1T Ψ−1 ψ(x))2 σy2 (x) = ρ2 σc2 + σd2 − ψ(x)T Ψ−1 ψ(x) + , (3.28) ψ(x)T Ψ−1 ψ(x) where µ = 1T Ψ−1 yT (1Ψ−1 1)−1 , y = {ycT , yeT }T . 3.4 Demonstrations of multi-fidelity surrogate The efficiency and effectiveness of the proposed MFS framework is illustrated by approxi- mating vascular G&R simulations of AAA expansion. Specifically, we first provide a demo study using univariate vascular simulations of AAA growth, aiming to show the improve- ments achieved by cokriging-based MFS compared to Kringing-based surrogate. Next, we demonstrate the practical use of the proposed vascular MFS using multi-fidelity G&R com- putational models. In this study, we employ a 3D G&R model in Section 2.3 as the HFM and an axisymmet- ric G&R model in Section 2.4 as the LFM. Dirichlet boundary conditions of both G&R models are prescribed by fixing the displacement of lower and upper boundaries of arterial geometries. A comparison between the two models are provided in Table. 3.1. Note that the computation time of all simulations are recorded on the same workstation with 10 cores i7-950 @3.3 GHz, 64GB RAM using the same sets of parameters. In addition, the time of the simulated stress-mediated adaptation process are fixed as 5300 days for all simulations. The cokriging method is implemented using a MATLAB toolbox ooDACE [Couckuyt et al., 2014]. 40 Table 3.1: Comparison between two-fidelity G&R models. 3D membrane G&R model (HFM) Axisymmetric G&R model (LFM) Software FEniCS MATLAB # finite elements 2528 25 Element type Linear Lagrange triangular element 1D quadric element Dimension 3D 2D (axisymmetry) Computation time a 70-90 minutes 10-15 minutes a The computation time denotes the estimated computational time of a G&R simulation 3.4.1 A univariate demonstration of MFS To demonstrate the advantages of the cokriging-based MFS compared to the kriging-based surrogate, we present an univariate demo of the AAA expansion. In particular, we take the length along the AAA centerline s as the only input variable for both HFM (3D membrane G&R) and LFM (axisymmetric G&R). All the other parameters are fixed by constant values so that we could focus on studying the simulations output, i.e., maximum diameters along the AAA centerline d, versus the input s. The simulated time is prescribed by t = 3000 days and the simulations are performed on an idealized geometry. Note that this demo is for illustration only and is not directly related to the rest of this section. In this demonstration, we randomly pick up a set of 6 HF sample points that is directly used to train a kriging-based surrogate. As Fig. 3.1 shown, the kriging, which is represented by the green dotted line, provides a smooth Gaussian shape approximation but is associated with large error under the area of s where less samples are provided. In contrast, the cokriging-based MFS is trained by 6 HF samples points as well as 50 cheap LF samples points. The MFS approximation is represented by the dark solid line in Fig 3.1, which shows a significantly improvement over kriging-based surrogate. This is because the GP term Zd () in cokriging formula improves the approximation accuracy in these regions with less HF samples. 41 Figure 3.1: A simple demonstration to cokriging-based MFS. s and d denote the length along the centerline and the associated diameter, respectively. The red and blue dashed lines represent the physics-based simulation outputs from HFM and LFM, respectively. The kriging-based surrogate is trained from 6 HF samples while the cokriging-based MFS is trained from the same set of 6 HF samples as well as 50 additional LF samples. 3.4.2 A MFS-based approximation to the G&R models of AAA growth As illustrated before, we can collect HF and LF data that incorporate the target values d under sets of four parameters, i.e., x = {θT , ξsT }T = {Kg , σdmg , s, t}T . Here, x can be divided into two subsets including a patient-specific intrinsic parameter set θ = {Kg , σdmg }T and a spatial-temporal parameter set ξs = {s, t}T . Essentially, θ includes the input parameters of the physics-based simulations, while ξs can be considered as the spatial-temporal locations in the output physical fields such as stress field and aortic diameters along centerline. Ac- cordingly, a single vascular simulation can generate multiple samples with the same value of θ but different values of ξs . Given the different properties of θ and ξs , a two-step sampling approach is developed to 42 70 2 1 0.5 4 1.2 70 60 50 40 1.2 4 0.51 2 3.5 65 3 1 60 1 2.5 dmg 55 dmg 2 0.8 50 0.8 1.5 45 0.1 0.6 0.6 0.2 1 40 65 55 45 0.2 0.5 4 35 0.065 0.07 0.075 0.08 0.065 0.07 0.075 0.08 Kg Kg (a) A contour plot of the approximated maxi- (b) A contour plot of the uncertainties esti- mum diameters. mated by MFS. Figure 3.2: Contour plots of the approximated maximum diameter and the associated un- certainty over the map of Kg and σdmg at t = 4000 days. The red dots represent the HF samples. The unit of color bar is mm. generate training data for MFS efficiently. First, we use a stratified sampling method to randomly sample 20 values of θ, and then input these values of θ independently into the HFM to collect HF simulation outputs. Second, we further randomly select 15 values of ξs from each HF simulation output. overall, we collect a set of Ne = 300 HF samples to train MFS. Meanwhile, the same approach is applied to the LFM to obtain a set of Nc = 1000 LF samples. HF and LF samples are used together to train the MFS. Fig. 3.2a shows a contour map of approximated maximal diameter under different values of patient-specific parameters at t = 4000 days, where the maximum diameter is defined as the maximal value of the diam- eters with all possible values of s. Additionally, Fig. 3.2b provides the standard deviation associated with the approximated maximal diameters. It shows that the uncertainty has been successful minimized under regions with sufficient training data. In these plots, the HF samples are represented by red dots, and each point on the map contains a set of 15 samples with different values of ξs . 43 80 MFSM-based approximation 80 MFSM-based approximation 95% confidence interval 95% confidence interval 60 60 d (mm) d (mm) 40 40 20 20 1000 2000 3000 4000 5000 0 20 40 60 80 t (days) s (mm) (a) Diameter against time (b) Diameter against length 80 80 d (mm) d (mm) 60 60 40 40 MFSM-based approximation MFSM-based approximation 95% confidence interval 95% confidence interval 20 20 0.065 0.07 0.075 0.08 0.4 0.6 0.8 1 1.2 Kg dmg (c) Diameter against Kg (d) Diameter against σdmg Figure 3.3: The MFS approximations of diameters against independent input variables in which irrelevant input parameters are fixed. The approximated error is indicated in blue shades. Now, given arbitrary values of input parameters x = {Kg , σdmg , s, t}T , the MFS is able to predict the target output, i.e., dˆ = ŷ = fsm ({Kg , σdmg }T , {s, t}T ), within 0.001 second. The collection of all target outputs is denoted as in silico data. Fig. 3.3 includes the four univari- ate plots that visually illustrate the relationship between the MFS-based approximations of diameters and the four inputted parameter. The uncertainty associated with MFS approx- imation is plotted by shaded confidence intervals. Note that the unrelated parameters are fixed in each demonstration. In the Fig. 3.3c and Fig. 3.3d, we notice that Kg strongly affects diameter d while σdmg 44 has less effect on diameter d. This is because Kg directly changes the growth rate, while σdmg mainly affects the shape of aneurysm. Furthermore, due to the low diameter-related uncertainty (the 95% confidence interval is less than 5mm), we conclude that the cokriging- based MFS is able to approximate the simulation with relatively small errors. 3.5 Discussion The PCM method significantly improves the efficiency of data generation from physics mod- els. In this study, it generates tens of thousands of IMDCs in 20 seconds, while computing a LF vascular G&R simulation takes 30 minutes. This means that it would take more than a year to provide the same amount of training data directly using a vascular G&R model. In conclusion, PCM is a powerful data augmentation tool which provides a basis for training data-driven models such as deep learning. The cokriging-based MFS is developed to generate an efficient approximation to a HF G&R model. The MFS can take advantages of multi-fidelity G&R models. In particular, a massive amount of computational cheap data simulated by LF vascular model merits the approxi- mation by greatly reducing the computational time of generating HF data. The success of the proposed MFS provides a solid foundation for more complex and time-consuming applications such as parameter estimation, validation, model reduction and patient-specific prediction, which will be discussed in later Chapters (Chapter 4, Chapter 5 and Chapter 6). Additionally, since the proposed MFS is a general simulation tool, it holds promise to approximate more types of cardiovascular disease simulations. 45 Chapter 4 AN AUTOMATIC PATIENT-SPECIFIC PARAMETER ESTIMATION FRAMEWORK USING MULTI-FIDELITY SURROGATE 4.1 Introduction The proposed MFS dramatically increases the efficiency of forward vascular simulations. One of the most valuable applications of MFS is to enable efficient and accurate patient- specific parameter estimation. Parameter estimation is generally considered to be an inverse problem, in which the patient-specific parameters should be carefully selected to minimize the discrepancy between simulation outputs and clinical observations. Nevertheless, this is 46 hard task that requires a novel MFS-based optimization framework that can automatically and efficiently solve the inverse problem of parameter estimation. Fortunately, optimization techniques have been broadly used in vascular modelings [Abraham et al., 2005, Marsden, 2013, Marsden et al., 2004, 2007, Quarteroni and Rozza, 2011]. In particular, the opti- mization approaches can be divided into two groups including gradient-based approach and derivative-free approach. In this study, we use the MFM, which cannot provide derivative information, to provide the basis for the parameter estimation, so we need to develop an optimization framework based the derivative-free evaluations. Recent studies have shown the promise for the derivative-free optimization using cokriging. Forrester et al. [2007] combined cokriging with an additional optimization algorithm which improves the accuracy of overall approximation by maximizing the improvements by re- samplings. Perdikaris et al. [2015] proposed a recursive cokriging-based framework to perform efficient design optimization under uncertainty. Additionally, a Bayesian optimization has been proposed to estimate resistance parameters in a artery bifurcation model based on cokriging [Perdikaris and Karniadakis, 2016]. Those studies provide illustrative results but cannot be directly applied for vascular parameter estimation using patient-specific clinical data. Therefore, inspired by those studies, we propose a new cokriging-based optimization framework particularly suitable for patient-specific parameter estimation of vascular models. However, there are two challenges remain to be addressed. First, the in silico and clinical data are required to be matched and compared following an automatic algorithm to enable the automatic parameter estimation framework. Second, approximation comes with uncer- tainty (3.28), which reduces the accuracy of parameter estimation. To cope with this issue, we propose an innovative algorithm to effectively search for the patient-specific parameters and reduce the related uncertainty in an iterative manner. In step 1, the MFS-based data, which denote the MFS approximations of G&R simulations, are assimilated with clinical data, which denote the patient-specific data measured in clinic. In step 2, an algorithm 47 of uncertainty quantification is developed to estimate the uncertainty associated with the discrepancy between the MFS-based data and the clinical data. In step 3, a Bayesian op- timization approach, i.e., Lower Confidence Bound (LCB) approach, is utilized to optimize the locations of sampling by balancing the exploitation of the extreme value of discrepancy and the exploration of the high uncertainty. Here, the LCB method is inspired by an optimization algorithm that treats the task of optimization as a multi-armed bandit problem [Srinivas et al., 2012]. Specifically, in Choi et al. [2008], a searching strategy switching between the exploration, i.e, searching for the sample location of with the largest uncertainty, and exploitation, i.e, searching for the sample location of current extreme, was designed for global maximization problem of swarming robot. Kahn et al. [2015] proposed a new sampling method, confidence bounding function, for swarm to reach a trade-off between the exploration and exploitation. In Srinivas et al. [2012], a intuitive Gaussian process upper confidence bound (GP-UCB) method was further performed, which shown a superior convergence rate for global optimization. In this study, inspired by these methods, the LCB approach is designed and implemented for solving our global optimization problem. The overall schematic flow of the iterative parameter estimation method is shown in the Fig. 4.1, in which the MFS is created in Section 3.4.2. The detailed algorithm is illustrated in the Algorithm 1. The method is implemented in MATLAB (Mathwork, Natick, USA). 48 Low-fidelity High-fidelity data data Construct MFS MFS-based Clinical data data Run new Step 1: high-fidelity Assimilated Data Assimilation simulation MFS-based and add the data output to Step 2: high-fidelity Uncertainty Quantification data Optimization objective with its uncertainty Step 3: Optimal Sampling with Exploitation and Exploration Optimal Sampling point Converge Figure 4.1: Schematic flow of the proposed MFS-based parameter estimation algorithm. 49 Algorithm 2 MFS-based parameter estimation for arterial G&R models 1: Get clinical data. 2: while the convegence is not reached do 3: Train MFS using HF and LF data (the simulation outputs from HFM and LFM, respectively). 4: Get MFS-based data by running MFS. 5: for ∀θ ∈ Ωθ do 6: Match MFS-based data and clinical data by using data assimilation algorithm. 7: Compute error associated with the MFS approximation. 8: Compute objective function J(θ) ˆ and related uncertainty σJ2 (θ) by using Monte- Carlo approach. 9: end for 10: Compute LCB surface Clcb (θ) from J(θ) ˆ and σ 2 (θ). J 11: Find optimal sampling location θ∗ through minimizing Clcb (θ). 12: Run HFM at θ∗ to get a new HF sample, and add it to HF data set. 13: end while 14: Let the optimal estimation of parameters θ opt = θ ∗ at the last iteration. 4.2 Patient-specific parameter estimation framework 4.2.1 Data assimilation To facilitate direct matching between the simulation outputs and the clinical data, we design a data assimilation approach that initializes spatial-temporal conditions [Balocco et al., 2010, Sermesant et al., 2006]. Consider a surrogate with two sets of input parameters x = {θT , ξ T }T , where θ denotes a Nθ dimensional vector of model parameters that controls the mechanical and geometrical features of a blood vessel, and ξ denotes a Nξ dimensional state variable, such as spatial location s and time t. The state variable in simulations and clinical measurements are defined in different coordinate systems. Accordingly, we specify the state variables of MFS-based data and clinical data by different notations, i.e., ξs and ξp , respectively. Given any model parameter θ, the MFS-based data can be represented by ŷ = fsm (θ, ξs ) that are simplified from Eq. (3.27). Meanwhile, the clinical data are also expressed in a similar form by using function y = fp (ξp ), illustrating 50 that a scalar clinical measurement y, such as arterial diameter, is measured at certain spatial- temporal locations represented by the state variable ξp . (1) (Ns ) A set of Ns clinical data consists of state variable Ξp = {ξp , ..., ξp } and the associated (1) (Ns ) measurements {fp (ξp ), ..., fp (ξp )}. For example, a set of clinical data may include a set of arterial diameters measured at different times and positions for a patient. In addition, a (1) (Ns ) set of MFS-based data also consists of state variables Ξs = {ξs , ..., ξs } and the associated (1) (Ns ) simulation outputs {fsm (θ, ξs ), ..., fsm (θ, ξs )}. Although the state variables Ξp and Ξs are measured in different coordinate systems, their elements should have the same interval values. For example, given two clinical measurements (1) (2) with state variable ξp = 1 cm and ξp = 2 cm, we can compute the interval ∆ξ (1) = 1 (2) (1) cm, and specify the interval value of the state variables of MFS-based data as ξs − ξs = ∆ξ (1) = 1 cm. Knowing the relative intervals of Ξs , we can anchor the value of an element in Ξs to complete data assimilation. Particularly, we choose to estimate the first state (1) variable ξs out of the whole set of MFS-based data Ξs by comparing the MFS-based data and clinical data. As a result, the data assimilation problem can be regarded as a simple (1) optimization problem, in which the ξs is estimated by minimizing this discrepancy between the MFS-based data and the clinical data, i.e., v u u 1 X Ns  2 (i) (i) ξˆs (θ) = arg min t (1) fp (ξp ) − fsm (θ, ξs ) , (4.1) (1) ξ ∈Ω Ns i=1 s ξ where Ωξ represents the constrained parameter space of ξ. Accordingly, the assimilated MFS-based data represented by ξˆs under any given θ is provided by ξˆs (θ) = ξˆs(1) (θ) + ξp − ξp(1) . (4.2) 51 4.2.2 Objective function and uncertainty quantification The assimilated state variable of the MFS-based data ξˆs (θ) is formulated by Eq. (4.2). Now, we focus on estimating the model parameters θ through an optimization approach. The objective function J(θ) of the optimization is quantified by v u u 1 X Ns  (i)  (i)  2 J(θ) = t fp ξp − fsm θ, ξˆs (θ) , (4.3) Ns i=1 which represents the discrepancy between assimilated MFS-based data and clinical data under the only independent variable θ. Eq. (3.28) shows that MFS-based data have approximate error that brings uncertainty to objective function J(θ). A Monte-Carlo method is used to quantify this uncertainty. Note that we use index α to distinguish different Monte-Carlo samples while using index i to represent different samples in clinical data set and MFS-based data set. Now, for any given model parameter θ with assimilated state variable ξˆs (θ), we can randomly sample a set of Monte-Carlo samples from the Gaussian distribution with the mean value ŷ({θT , ξˆsT (θ)}T ) and the variance σy2 ({θT , ξˆsT (θ)}T ) provided by MFS referring the Eq. (3.27) and Eq. (3.28). (α) (i) Then we can represent the αth Monte-Carlo sample in a function form, i.e., fˆsm θ, ξˆs (θ) ,  where α = 1, ..., Nm . The value of the objective function that corresponds to the αth Monte- Carlo sample is then provided by v u u 1 X Ns  (i)  (α) (i) 2 (α) J (θ) = t ˆ ˆ fp ξp − fsm θ, ξs (θ) . (4.4) Ns i=1 By collecting all J (α) (θ) of Monte-Carlo samples, the objective function and the associated 52 uncertainty can be approximated as Nm  ˆ 1 X (α)  J(θ) = J (θ) , (4.5) Nm α=1 Nm  1 X ˆ 2 σJ2 (θ) = (α) J (θ) − J(θ) , (4.6) Nm − 1 α=1 respectively. In practice, we divide the domain Ωθ into grids, and then estimate the objective functions and the associated uncertainty at the grid points. Although this algorithm requires a large number of forward simulations, the overall computational efficiency is still high due to the use of MFS. In particular, the cokriging-based MFS is capable of generating a random field by performing 1000 approximations at different values of θ within 0.2 seconds in a 3.3 GHz CPU core. 4.2.3 Lower confidence bound-based optimal sampling The estimated discrepancy formulated in Eq. (4.5) provides a solid foundation for optimal parameter sampling. Nevertheless, during the parameter searching, the uncertainty in the parameter space shown in Eq. (4.6) can be large and a direct minimization of discrepancy could lead to large error. We therefore need a sampling method that iteratively update the MFS by adding new samples to find the best estimated parameters. To achieve this goal, we adopted a combined strategy that balances both exploitation and exploration. In particular, we seek a new measurement that samples point by trading off between exploiting the local minima and exploring the high-uncertainty in the global field. This strategy is based on a Bayesian optimization approach, i.e., Gaussian Process Upper Confidence Bound (GP-UCB), in which an intuitive upper bound of the possible maximum over the field was developed [Choi et al., 2008, Kahn et al., 2015, Srinivas et al., 2012]. In 53 this study, we replace the upper bound with a lower bound denoted by LCB surface Clcb (θ), and the optimal sampling location is given by ˆ − p θ∗ = arg min Clcb (θ) = arg min J(θ) βσJ (θ), (4.7) θ∈Ωθ θ∈Ωθ where the coefficient β is chosen according to the structure of the problem to ensure the convergence. The LCB surface Clcb provides a reasonable lower bound that balances both ˆ effects of small mean value J(θ) and large uncertainty σJ (θ). Srinivas et al. [2012] showed that this algorithm is promised to be globally converged if it operates on a function with a high level of smoothness. Up to this point, we have formulated the framework from the MFS training to the optimal sampling. Next, the MFS should be iteratively trained using the augmented HF sample data Xe , which contains both original HF data set and the additional data samples found by the LCB sampling algorithm. The whole process is repeated until the convergence criterion is met. The optimal fitted patient-specific parameter is denoted by θopt . 4.3 Case study setup and clinical data In Section 4.3.1, we explain the details of the clinical data of AAAs. In Section 4.3.2, the G&R and MFS are set up to capture AAA expansions efficiently. 4.3.1 Clinical data of AAAs 82 follow-up CT images of AAAs from 21 patients are retrospectively obtained at Seoul Na- tional University Hospital. Table. 4.1 summarizes the demographic information of patients, the number of follow-up scans, the time of follow-up scans, and the maximum diameters of 54 AAAs during each scan. The average and median age of patients are 66.9 and 67 years old, respectively. The average number of CT scans is 3.9. Table 4.1: Demographic data and CT scan information of individual patients. Patient ID # Scans Gender Age Time of scans (days) Max. diameters (mm)a P01 3 Male 71 0, 203, 733 42.9, 44.8, 51.1 P02 2 Male 71 0, 632 37, 39.1 P03 3 Female 64 0, 494, 1357 47.9, 54.4, 65 P04 5 Male 65 0, 182, 361, 538, 728 41.7, 42.8, 45.7, 48.7, 51.1 P05 5 Male 74 0, 347, 702, 1054, 1223 48.1, 50.2, 52.2, 57, 58.4 P06 5 Male 66 0, 374, 1074, 1438, 2136 30, 31.7, 33.2, 33.9, 34.8 P07 5 Male 54 0, 386, 757, 1121, 1290 38.9, 40.6, 42.5, 44.9, 46.8 P08 5 Male 62 0, 227, 674, 1049, 1403 38.9, 41.3, 43.4, 46.8, 54.5 P09 3 Male 73 0, 97, 266 40.5, 41.6, 44.4 P10 4 Male 59 0, 522, 922, 1344 39.6, 41.5, 43.4, 46.9 P11 4 Male 54 0, 399, 774, 1152 40, 43.4, 46.2, 50.6 P12 3 Male 78 0, 523, 873 41.7, 44.8, 46.6 P13 4 Male 68 0, 543, 691, 874 39.2, 44.9, 46.9, 48.9 P14 3 Male 71 0, 349, 714 43, 45.2, 51.5 P15 5 Male 67 0, 183, 366, 534, 709 44.9, 46.2, 47.7, 49.7, 52.3 P16 3 Male 72 0, 189, 526 45.4, 46.7, 48.2 P17 5 Male 72 0, 246, 421, 587, 783 37.2, 41.8, 41.8, 45.2, 47.4 P18 4 Female 65 0, 613, 1048, 1515 32.3, 39.2, 45, 50.7 P19 4 Male 78 0, 309, 1976, 2310 34.7, 36.9, 54.2, 60.8 P20 4 Male 64 0, 729, 922, 2359 36.4, 40.4, 41.6, 50.9 P21 3 Male 57 0, 440, 999 38.7, 43.2, 54.2 a Max. diameters includes the maximum diameters of AAAs during each scan. To compare the morphological data from patients and the MFS simulation results, we utilize an Inscribed Maximal Diameter Curve (IMDC) Method [Gharahi et al., 2015] that can compute the aortic diameter along its centerline. The IMDC approach first selects the part of an AAA image ranging from the iliac bifurcation to the renal branches, and then measures the inscribed sphere diameters by moving the inscribed maximal sphere through the aorta from the bottom to the top. Accordingly, a curve of maximum diameter (d) against the axial location along the arterial centerline (s) can be generated for every CT scan image using the IMDC method. In addition, each patient has multiple CT scans captured at different times from a follow-up study, so each set of morphological data is labeled at a given time point (t). 55 As a result, the morphological data of patients can be transformed and rearranged into a function d = fp (ξp ) that represents the diameters measured at the state variable ξp = {s, t}T . This is consistent with the corresponding content shown in Section 4.2. 4.3.2 G&R models and MFS setup Computational vascular G&R models (Chapter 2) are employed to simulate the morpho- logical changes of AAAs by capturing elastin degradation and collagen adaptation in the aortic wall. Elastin is a main load-bearing elastic constituent in arterial tissue, but worn or damaged structural elastin cannot be reproduced in adulthood. A large portion of elastin degradation causes an increase in local stress that may lead to dilation of the aorta [Wilson et al., 2012]. Also, Rizzo et al. [1989] found that elastin degradation is a key feature that is ubiquitous in AAAs. Therefore, we choose the parameter σdmg , which determines the shape of elastin degradation function, as one of the model parameters to be estimated. More details of elastin degradation are introduced in Section 2.2.4. In the constrained mixture model, the parameter Kg , which affects the turnover rate of collagen fibers, has a major impact on the growth rate of AAA expansion (Section 2.2.3). Thus, Kg is selected as the second model parameter to be estimated. In summary, we employ two key model parameters, i.e., θ = [θ1 , θ2 ] = [Kg , σdmg ], to simulate AAA enlargement, whereas the other parameters in G&R models are fixed by population-based values, as given in Seyedsalehi et al. [2015], Zeinali-Davarani et al. [2011]. Now, given four parameters including two model parameters and two state variables x = {θT , ξsT }T = {Kg , σdmg , s, t}T , we can obtain two values of d from HFM and LFM, respec- tively. As aforementioned in Section 3.4, a FEniCS based 3D G&R model and a MATLAB based 2D axisymmetric G&R model introduced in Section 2 are employed as the HFM and the LFM, respectively. The same IMDC technique presented in Section 4.3.1 is used to con- 56 vert the 3D G&R simulation outputs into morphological curves that can be directly matched with the axisymmetric G&R simulation outputs. 4.4 Results We use a two-step sampling approach to obtain a set of Ne = 300 HF samples a set of Nc = 1000 LF samples. The MFS is then trained by using both HF and LF samples. The details of MFS results are described in Section 3.4.2. In summary, given arbitrary values of input parameters x = {Kg , σdmg , s, t}T , the MFS is able to predict the aortic diameter, i.e., dˆ = fsm ({Kg , σdmg }T , {s, t}T ), within 0.001 second. The clinical data of AAAs are introduced in Section 4.3.1. The patient-specific model pa- rameters of G&R models are estimated by using the optimization framework shown in Sec- tion 4.2. Fig. 4.2 provides an example (Patient 18) of the iterative sampling, in which the LCB surface Clcb (θ) at the final iteration is plotted over a local domain of θ. Each red dot represents a set of newly added HF samples that are used to update MFS. The black arrows indicate the order of the added sample points, illustrating how the LCB approach greedily searches for the best-fit parameters. The LCB approach finally converges to the region with low values of Jˆ near the upper left corner. Meanwhile, it also explored the points near the lower right corner to minimize large uncertainties. The re-sampling and updating of MFS have been repeated 10 times prior to the convergence. A MFS-based personalized G&R simulation of AAA is obtained by imposing the best-fitted patient-specific parameters into MFS. Fig. 4.3 shows the AAA curves of Patient 18 obtained from the individual G&R simulation and clinical data. Red dashed lines with circles repre- sent the clinical data with different values of state variable. Black curved lines represent the patient-specific simulation outputs. Blue shaded regions represent the estimated approxi- 57 Figure 4.2: A contour plot that provides the LCB surface Clcb (θ) at the final iteration of optimization. The unit of color bar is mm. The red dots denote the additional samples used to iteratively update the MFS. The black arrows indicate the order of additional samples. 58 60 Patient-specific clinical data Patient-specific simulation Uncertainty 50 t = 4.15 years t = 2.87 years d (mm) 40 30 t = 1.68 years t = 0 year 20 20 30 40 50 60 s (mm) Figure 4.3: Diameter (d) versus length along centerline (s) of the AAA of Patient 18. The red dashed lines with circles represent the clinical data measured from follow-up CT images captured at t = [0, 1.68, 2.87, 4.15] years. The black curves represent MFS-based patient- specific G&R simulation outputs. The uncertainty related to MFS approximation is indicated by 95% confidence interval and is shown by blue shades. mation error between the MFS and original high-fidelity G&R. The region of approximation error is, however, too small to be visible in the figure. The entire process of parameter estimation is repetitively conducted for each of the 21 pa- tients. Table. 4.2 records the estimated parameters, computational time and number of iterations used in the optimization for each patient. Note that the overall computational time includes the time spent on MFS model training, data assimilation, the generation of random fields and the LCB searches. The generation of HF and LF data has been done all ˆ opt ), which at once before parameter estimation. In addition, the final objective function J(θ represents the discrepancy between the personalized simulation outputs and clinical data, is 59 small (less that 2 mm) compared to AAA diameters (ranges from 20 mm to 70 mm). This illustrates that the proposed patient-specific simulations are capable of fitting the follow-up images of AAAs. The associated uncertainty, i.e., σJ (θopt ), for each of the 21 patients is ˆ opt ), indicating that the proposed framework is successfully minimized to less than 5% of J(θ robust in estimating the patient-specific parameters among different patients. The compu- tational time for each parameter estimation, which ranges from 5.1 to 64.3 minutes, is less than the clinical-relevant time (1 -2 hours) [Taylor et al., 1999]. The computational time has no significant correlation with variables such as the number of CT scans and the growth rate of AAAs. Table 4.2: Information of patient-specific parameter optimization. Patient ID Kg σdmg J(θˆ opt ) σJ (θopt ) #I a tb P01 0.0696 0.9331 0.5589 0.0123 6 11.9 P02 0.0815 1.0364 0.4543 0.0168 15 33.5 P03 0.0688 0.9962 1.6192 0.0472 7 17.0 P04 0.0676 1.0974 0.6816 0.0221 2 4.0 P05 0.0699 1.1446 1.4582 0.0128 13 42.8 P06 0.0823 0.4154 1.1580 0.0503 7 23.6 P07 0.0774 0.7769 0.8491 0.0130 6 22.4 P08 0.0697 1.0846 1.4228 0.0095 19 64.3 P09 0.0605 0.9667 0.9239 0.0396 4 12.1 P10 0.0792 0.4513 1.0106 0.0440 4 11.6 P11 0.0743 0.9385 0.6057 0.0094 6 18.4 P12 0.0778 0.9744 0.5712 0.0134 12 37.5 P13 0.0701 1.2372 0.8478 0.0175 13 49.4 P14 0.0701 0.8654 1.2773 0.0101 5 16.5 P15 0.0699 1.1962 0.7903 0.0361 9 32.2 P16 0.0787 0.7744 0.5394 0.0224 12 36.6 P17 0.0693 0.9154 0.9867 0.0140 2 5.1 P18 0.0673 0.6359 0.9742 0.0301 10 28.0 P19 0.0697 1.0769 1.1565 0.0382 7 23.8 P20 0.0759 1.0154 1.1671 0.0370 6 21.5 P21 0.0694 0.4615 0.9649 0.0207 15 51.5 a b #I is number of iteration. t with unit of minutes is the computation time of parameter optimization. The population-based values of Kg and σdmg are collected from the 21 patients. The sample mean and the standard deviation of Kg is 0.0723 and 0.0055, respectively; and the sample 60 mean and the standard deviation of σdmg are 0.9045 and 0.0575, respectively. This informa- tion can be employed as prior distributions of parameters in future studies. 4.5 Discussion This study has made contributions to the vascular G&R modeling field by establishing a new physics-based machine learning framework, capsulized as (1) improving the computational efficiency of arterial G&R simulation, and (2) providing an efficient G&R-based optimiza- tion workflow to estimate personalized parameters from clinical data. The computational efficiency of the overall simulation framework is significantly improved without compromis- ing on accuracy and a good balance is achieved between fundamental laws of physics and computational efficiency. The new framework is demonstrated by using clinical CT images of AAAs from 21 patients. Significant contributions of this work are discussed here. The primary contribution of this work is to greatly improve the computational efficiency in G&R based simulations, hence providing fast estimations of patient-specific parameters and individual predictions of arterial disease progression. The significant improvement in computational efficiency is mainly attributed to a variety of novel techniques. First, the computational time of a high-fidelity simulation (a three dimensional finite element G&R simulation) is reduced from roughly one week [Zeinali-Davarani et al., 2011] to 70 -90 minutes by implementing the FEniCS [Alnæs et al., 2015]. Second, a large number of cheap simula- tions (axisymmetric finite element G&R model [Baek et al., 2006]) are conducted to remedy time-consuming high-fidelity simulations. Third, a cokriging-based MFS is constructed as an important intermediate stage connecting the forward arterial simulation to the inverse problem of personalized parameter estimation, and each simulation of the MFS only takes around 0.0002 seconds. Finally, the LCB optimal sampling approach based on the Bayesian optimization [Choi et al., 2008, Kahn et al., 2015, Srinivas et al., 2012] is developed to further 61 improve the efficiency of estimating the personalized values of intrinsic parameters so that each parameter estimation can be conducted within the clinically relevant time frame (1-2 hours) [Taylor et al., 1999]. Therefore, the above techniques enable us to effectively overcome the current challenge of the inefficiency in the computational G&R modeling [Ambrosi et al., 2019], thus laying a foundation for further expanding the applications of G&R models. Another major contribution of this work is to develop a novel iterative multisteps optimiza- tion algorithm that enables personalized parameter estimation of G&R models (illustrated in Fig. 4.1). The proposed optimization algorithm consists of three main steps, including data assimilation, uncertainty quantification, and a Bayesian optimization based sampling. By performing these three steps iteratively, the discrepancy between G&R simulation outputs and clinical data can be effectively minimized while avoiding local minima. Particularly, the first step allows us to assimilate G&R simulation outputs and clinical CT images following an objective standard, and the second step allows us to quantify the uncertainty associated with the discrepancy between the assimilated MFS data and clinical data via a Monte-Carlo approach. These two steps together provide an important basis for the third step, i.e., the optimal sampling via the LCB algorithm, that aims to balance the exploitation of extreme and exploration of the high uncertainty. This ensures a fast global convergence with a high degree of smoothness. The proposed optimization is similar to other kriging-based optimiza- tion methods such as Efficient Global Optimization (EGO) that has been broadly applied to find the global extreme by balancing the local and global searching [Forrester et al., 2007, Perdikaris and Karniadakis, 2016], as well as the surrogate management framework [Booker et al., 1998] that contains two distinct steps, i.e., search and poll, for a fast convergence. These two optimization methods, however, cannot be applied to this study, because they re- quire the optimization target to follow GP random field and only guarantee the convergence to a local minimum, respectively. As far as we know, this is the first study that systematically estimates the patient-specific 62 values of parameters Kg and σdmg from follow-up CT images. Here, Kg and σdmg are two significant intrinsic parameters that control the progression of AAAs. Although there exists a large number of studies on parameter estimation [Marsden et al., 2008, Sankaran et al., 2013, Szafron et al., 2019], most do not estimate personalized intrinsic parameters by using clinical data. In our previous studies [Jiang et al., 2020, Zhang et al., 2019], where these intrinsic parameters were adjusted to generate in silico data for training the predictive models, the limited knowledge of these intrinsic parameters confined the accuracy of the prediction. However, with the estimated values of Kg and σdmg that can take personalized variability into account, the simulation results will become more realistic and thus enhance predictability. The predicted diameters for 17 out of 20 patients fell within the tolerance interval. The average prediction error of the diameter and the maximum diameter along centerline are 1.39 mm and 1.67 mm, respectively, which are smaller than the prediction error obtained by previous studies [Akkoyun et al., 2020, Jiang et al., 2020, Lee et al., 2018]. In summary, the successful parameter estimation and prediction of AAA expansion demonstrate that the MFS holds promise for effective parameter estimations and predictive medicine. For the sake of computational efficiency, this study has to limit the number of parameters. In particular, we focus on estimating four key parameters which are acceptable in capturing the time-dependent morphological features of AAA growth. In the future, we plan to uti- lize high-dimensional optimization techniques, such as Sarkar et al. [2019], to include more parameters that can affect AAA growth, such as arterial tortuosity [Akkoyun et al., 2020] and hemodynamics-related parameters [Zambrano et al., 2015], thus ensuring higher inter- personal variability. Another limiting factor is that we only quantify the modeling error between in silico data and the clinical data during the model training while the clinical mea- surement error is not taken into account. This is because this study primarily focuses on the application of physics-based machine learning. The measurement error can be conveniently integrated with the clinical measurement error in the future works (e.g., using the previous work of Bayesian calibration [Zhang et al., 2019]). Moreover, a set of 21 follow-up data may 63 limit the prediction performance, we plan to collect more high-resolution follow-up data to further improve the prediction performance in the future. The proposed framework holds promise for being extended to more applications. This is not only because of its high efficiency, but also because it is a general framework consisting of multiple submodels that can be easily replaced to enable different applications. For example, the arterial G&R models in the proposed framework can be replaced by cardiac G&R models [Klepach et al., 2012] to enable the fast surrogate and personalized parameter estimation of a cardiac problem. Besides, the fidelity of the framework can also be further improved by using alternative physical models and machine learning models. For example, a shell element G&R model [Laubrie et al., 2019] can be used to improve the fidelity of MFS simulation. 64 Chapter 5 A PHYSICS-BASED DEEP LEARNING APPROACH TO PREDICT ABDOMINAL AORTIC ANEURYSM ENLARGEMENT 5.1 Introduction In recent years, notable advances in statistical tools have been implemented to predict the maximum diameter with longitudinal AAA scanning data. Sweeting and Thompson [2012] developed a hierarchical linear growth model utilizing a zero-mean Gaussian distributed random-effects term to simulate the growth effects of aneurysms. Others have used linear and quadratic hierarchical growth models to make predictions of the evolution of aneurysms [Brady et al., 2004, Eriksson et al., 2005]. Do et al. [2019] tested a method of dynamic Gaussian process to predict three-dimensional surface evolution and its uncertainty using 65 patient follow-up images. Nevertheless, due to the limited sample size and large measurement error of longitudinal data from patients, all these statistical predictive tools are not accurate enough to aid in clinical treatment. To cope with this problem, in this study, a novel predictive tool is designed using a deep learning algorithm, which holds promise for clinical treatment and recommendation. Deep learning and deep architectures in general have been applied in an enormous number of research areas, with the majority being in computer vision [Memisevic and Hinton, 2007] and natural language processing [Collobert and Weston, 2008]. Deep learning has also been widely applied in risk prediction based on electronic health records [Cheng et al., 2016], image labeling [He et al., 2004], traffic flow prediction [Huang et al., 2014], image segmentation [Long et al., 2015], medical image segmentation [Lai, 2015], stress estimation [Liang et al., 2018], and many other fields. Deep architecture has been investigated since 1980 [Fukushima, 1980] and proved to be more effective and requires fewer resources than a shallow structure of the same size, i.e., same number of nodes. The merits of deep structure come from its ability to improve efficiency by distributing different kinds of tasks through different layers [Lai, 2015]. For instance, the low layers can perform low level tasks like gradients computation or edge detection while the higher layers can perform classification or regression. However, as the networks are constructed in deeper layers, the training becomes prohibitively slow due to the problem of ‘vanishing gradients’ [Bengio et al., 1994]. In particular, when the error is back-propagated from the output layer, it is multiplied by the derivatives of the activation function, which is near zero for those saturation nodes. Consequently, the error, as the driving force for the gradient decent algorithm, is dramatically dissipated which results in an extremely slow training rate for those nodes behind the saturated node. The training problem remained until 2007 when Hinton proposed a two-stage learning scheme based on the Restricted Boltzmann Machine (RBM) [Hinton, 2007]. After that, other meth- ods, such as Rectifier Networks [Glorot et al., 2011], drop out technique[Hinton et al., 2012], 66 and Convolutional Neural Network [Meng et al., 2015] have been further developed and improved to solve the training problem. Research has shown that deep structure yields a high level of generalization and a low test error only when it is trained on a large training set. For instance, LeCun et al. [1998a] utilized a MNIST [LeCun et al., 1998b] dataset of hand-written numbers with 60, 000 samples to train a 3-layer deep network. Unfortunately, such a large dataset of longitudinal AAA images is unavailable. Therefore, the applications of deep structure in medical data are limited to medical image segmentation. The two main contributions of this work are as follows. First, to address the fundamental problem of the limited longitudinal data size, a massive in silico dataset is generated using a physics-based computational model and an approximation algorithm. Second, a novel predictive tool using a deep learning algorithm is established to combine both in silico and measured longitudinal data. To achieve these contributions, we employ a Deep Belief Network (DBN) to predict the AAA shape in a regression framework. To cope with the limited dataset, a small in silico dataset is generated by a computational Growth and Remodeling (G&R) model that can simulate the evolutionary process of AAAs. In the previous Section 3.2, the Probability Collocation Method (PCM) was introduced as an approximation method to reproduce a large amount of in silico data based on the G&R simulation outputs, which enables a computationally efficient data generation process. In Section 5.2, inspired by Hinton [Hinton, 2007], a two- stage learning scheme is employed to train the DBN. Briefly, the network is first pre-trained with the in silico data by the RBM in an unsupervised manner. The network is then fine- tuned with the labeled patient data. The data processing and model testing results are shown in Section 5.3, which is followed by the discussion and conclusion in Section 5.6. A basic schematic drawing of the overall flow is shown in Figure 5.1. To the best of our knowledge, this is the first effort to adapt a deep structure coupled with a G&R computational model to predict AAA enlargement. 67 Massive in silico G&R simulation PCM Two-stage outputs approximations model training data Small Follow-up CT Inscribed maximal patient Model testing images of AAAs diameter curves data Figure 5.1: Diagram of overall methodology. A massive number of in silico data and a small number of patient data are collected, which is followed by a two-stage model training and a model testing. 5.2 Deep neural network algorithm In this section, we introduce the constructing and training of DBN. A standard structure of the DBN Hinton et al. [2006] with two layers of RBM [Hinton, 2010] is utilized as shown in Figure 5.2. Assume that we have two types of variables: the visible unit (x) and the hidden unit (h). The two variables are governed by an energy function E(x, h). Given that both visible and hidden units follow binomial distribution, a Boltzmann Machine is defined as an energy-based model using a second-order polynomial [Bengio, 2009] E(x, h|θ) = −bT x − cT h − hT W x − xT U x − hT V h, where θ is the collection of b, c, W , U , and V . Thus, any probability density function (PDF) of visible layer P (x), as well as joint and conditional PDFs, can be easily represented by a 68 Figure 5.2: The deep architecture of the DBN. Two layers of RBM are trained in an unsu- pervised manner (pre-trained) using CD-1 algorithm. The top layer utilizes a neural network sigmoid regression for the prediction. normalized form of the energy function. For instance, the PDF of x can be computed as X e−E(x,h) P (x) = , (5.1) h Z e−E(x̃,h) is the normalization factor and x̃ can be all possible values of P P where Z = x̃ h the visible vector x. The realization of x̃ can be considered as reconstructed visible units. The parameter θ is estimated using the maximum likelihood estimation. However, the energy-based PDF in this study requires sampling of two conditional probabilities: P (h|x) and P (x, h). Although this can be done via the Markov chain Monte Carlo (MCMC) sam- pling [Hinton et al., 1984], it is highly computationally expensive. Therefore, contrastive divergence (CD), an alternative way of finding the log-likelihood, is utilized to provide a more efficient solution. 69 The RBM is introduced by posting an additional restriction: U = 0 and V = 0. In other words, there is no connection between the units in the same layer, either visible or hidden. Note that there is no link among units in the same layer in Figure 5.2. Furthermore, to incorporate the real-valued values, we specifically assume that the visible units follow Gaus- sian distribution, i.e., xi ∼ N (ai , σi ), and the hidden units follow binomial distribution, i.e., hj ∈ {0, 1}. Thus, the modified energy function is defined as [Hinton, 2010] V H V X H X (xi − ai )2 X X xi E(x, h|θ) = − − cj hj − wij hj (5.2) i=1 2σi2 j=1 i=1 j=1 σi The conditional PDF of the visible units given the hidden units can be computed as e−E(x,h|θ) P (x|h, θ) = P −E(x,h|θ) . xe Note that P (h|x, θ) can be computed in the same manner. Using E(x, h|θ) in (5.2), we have V ! X P (hj = 1|x, θ) = sigm wij xi + cj , i=1 H ! (5.3) X P (xi = x|h, θ) = N ai + σi hj wij , σi2 , j=1 1 where sigm(x) = 1+exp(−x) is the sigmoid function. The likelihood gradient can be computed by taking the derivative of P (x) in (5.1) with 70 respect to θ, and we have log P (x) ∂θ 1 X ∂E(x, h) = −P e−E(x,h) h e−E(x,h) h ∂θ 1 X X −E(x̃,h) ∂E(x̃, h) + e (5.4) Z x̃ h ∂θ X ∂E(x, h) X X ∂E(x̃, h) =− P (h|x) + P (x̃, h) h ∂θ x̃ h ∂θ     ∂E(x, h) ∂E(x̃, h) = −EP (h|x) + EP (x̃,h) . ∂θ ∂θ The expectation EP (h|x) [.] is also called positive phase distribution or data distribution while the other expectation EP (x̃,h) [.] is called negative phase or model distribution [Bengio, 2009]. Optimization of (5.4) involves sampling from P (x̃, h) and it can be realized by running a Gibbs sampling until it reaches equilibrium, which is extremely time consuming. Alterna- tively, Hinton [2002] suggested the CD learning which minimizes the difference between the data distribution and the one-step reconstructed distribution rather than directly minimizing the difference between the data and the model distribution. Empirical studies have shown that the CD method is efficient and effective enough to make the DBN unsupervised learning practical. Given the approximation to the derivative of P (x), we can pre-train the DBN by updating the weights iteratively using the in silico data in an unsupervised manner. The training is performed throughout the DBN structure shown in Figure 5.2. We consider this step, i.e., the pre-training of the DBN, as the first step of the two-stage learning scheme, enabling the DBN to capture the changes of geometrical features simulated by the G&R of AAA expansion. After the pre-training, the DBN is unfolded into a Neural Network (NN), which is further 71 trained in a supervised manner. This supervised learning, i.e., fine-tuning, is considered the second step of the two-stage learning scheme. Specifically, during the fine-tuning, the pre- trained weights of the unfolded DBN are properly adjusted for a better ability in capturing patient-specific features of aortic enlargement from the patient data. Our proposed two-step training model could be interpreted as a deep learning version of the Bayesian approach, where computer-generated data act as a prior distribution and the patient data for fine-tuning can be viewed as new measurements to compute the posteriori distribution for prediction [Lee et al., 2017a, Neal, 1996, Williams, 1996]. We also would like to remind that the Bayesian approach is normally used for prediction with the limited sized data (by leveraging the prior distribution), which is well suited for our case. 5.3 Data processing and results In this section, we introduce the data processing step and demonstrate the effectiveness of our proposed predictive model using observations of patient-specific CT images. 5.3.1 Data processing Given CT images of AAAs taken from a patient, we can obtain IMDCs with regular time intervals. Let ft,i be an IMDC of the ith AAA obtained at certain scan time noted as t. A collection of f generated at different times, i.e, t − 2, t − 1, t, and t + 1, provides us a timeline growth of the AAA. Note that the time intervals between IMDCs are fixed as the same constant. For each data point, we choose the feature vector xi as the collection of the three adjacent IMDCs, and the prediction target yi as the IMDC at the next time step in 72 the future, xi = [ft−2,i , ft−1,i , ft,i ] , yi = ft+1,i . Following Table 5.1, however, we realize that the time intervals between patient follow-up CT scans are not constants. To address this problem, a fixed-time interval of 270 days is chosen and all patient data are linearly interpolated by this fixed-time interval. Generally, we obtain 55 sets of interpolated patient data as collections of {xi , yi } from 20 patients. Specifically, 6 sets of the patient data from 6 different patients are randomly selected as testing data while the others are employed for pre-training. Table 5.1: Demographic data of patients Patient ID # Scans Gender Age Time of scans (days) P01 3 Male 71 [0, 203, 733] P02 3 Female 64 [0, 494, 1357] P03 5 Male 65 [0, 182, 361, 538, 728] P04 5 Male 74 [0, 347, 702, 1054, 1223] P05 5 Male 66 [0, 374, 1074, 1438, 2136] P06 5 Male 54 [0, 386, 757, 1121, 1290] P07 5 Male 62 [0, 227, 674, 1049, 1403] P08 3 Male 73 [0, 97, 564] P09 4 Male 59 [0, 522, 922, 1344] P10 4 Male 54 [0, 399, 774, 1152] P11 3 Male 78 [0, 523, 873] P12 4 Male 68 [0, 543, 691, 874] P13 3 Male 71 [0, 349, 714] P14 5 Male 67 [0, 183, 366, 534, 709] P15 3 Male 72 [0, 189, 526] P16 5 Male 72 [0, 246, 421, 587, 783] P17 4 Female 65 [0, 613, 1048, 1515] P18 4 Male 78 [0, 309, 1976, 2310] P19 4 Male 64 [0, 729, 922, 2359] P20 3 Male 57 [0, 440, 999] Given the multiple sets of G&R input parameters γ = [Kg , σd , µd ], the G&R model associated with the PCM approximation can produce a large number of longitudinal 2D profile curves 73 to capture the enlargement of AAAs during a time span, which can be transformed into the in silico data, i.e., artificially generated collections of {xi , yi }. In this study, we focus on predicting the aneurysm growth; hence, only those in silico data with maximum diameters ranging from 3 to 8.5 cm are accepted into the training dataset. After random samplings and rejections, 32900 sets of in silico data {xi , yi } are collected in the training dataset. In order to assimilate patient data and in silicon data together with the same dimensionality of data, we trim all both patient data and in silicon data into regions where the coordinate along the centerline ranges 0 to 8 cm. Next, each shape curve is discretized by a grid size of 81, meaning that the dimension of xi and yi are 243 and 81, respectively. Additionally, it has been shown that it is much simpler to train the RBM by the data with a zero mean and unit variance. [Hinton, 2010]. Thus, we normalize the training data before training the deep learning model. Specifically, a scale factor of 8.5 cm is selected for the normalization, based on the fact that the largest diameters of all our patient data are significantly smaller than 8.5 cm (also immediately recommended for surgical options [Brady et al., 2004]). As a result, all diameters obtained from both the patient data and the simulated results are normalized into the range [0, 1] before the training. As described in Section 5.2, the 32900 sets of normalized in silico data are utilized to pre- train the DBN in an unsupervised manner. Next, the selected patient data are employed to further update the deep structure in a supervised manner. Finally, during the model testing, we can collect the predicted IMDCs, which are then transformed back to the normal scale as the final prediction results. The overall method is depicted in Figure 5.3. 74 PCM outputs Real patient images In silico data 𝒙𝒊 In silico data 𝒚𝒊 Visible layer Hidden layer 1 DBN Hidden layer 2 Neural network Fine-tuned network Real patient data {𝒙𝒊 , 𝒚𝒊 } Make prediction Figure 5.3: Diagram of the model training. The DBN is pre-trained by in silico data in an unsupervised manner, and is unfolded into a neural network. Next, the neural network is fine-tuned with the patient data in a supervised manner to predict the AAA expansion. 5.3.2 Test set-up For a deep structure, parameters, such as the number of hidden units and the number of epochs1 , are important factors in determining the model performance. To avoid over-fitting, as a rule of thumb for the generative models using the high-dimensional data, the number of parameters is constrained [Hinton, 2010]. In our case, the data dimensionality (243) is 1 Epochs are the number of times that the model is trained through the whole training set. 75 significantly lower than the size of training samples (32900). In order to test the effect of the number of nodes within each layer, we construct a number of 2-layer DBNs with different sets of hidden nodes in each layer. In particular, since there is a remarkable difference in the dimensions of the data and the label, i.e., 243 versus 81, we test the DBNs of three structures: increasing width, decreasing width and equal width. As shown in Table 5.2, as the number of input or output layer nodes decreases, it becomes harder for the model to capture the representative features in the data. Note that in Table 5.2, the prediction error is quantified by the discrepancy between the predicted IMDC and the patient IMDC using the standard root mean squared error (RMSE), and the average RMSE is calculated from of 6 test samples. We then conclude from the test set-up that 300 nodes in each layer leads to the smallest prediction error among all tested configurations. Table 5.2: Effect of the number of nodes in a 2-layer DBN on the model testing. Number of nodes Training time Average RBM-1 RBM-2 (seconds) RMSE (cm) 1000 50 31 0.264 500 100 21 0.186 300 300 24 0.180 100 500 18 0.192 50 1000 23 0.2 As aforementioned, one of the problems in the DBN approach is that the pre-training gen- erates a large in silico dataset (32900 samples), while the patient data for fine-tuning are limited (49 samples). To better employ both datasets, we fix the epoch of the pre-training process to be 1 while changing the number of epochs of the fine-tuning process. The RMSE and training time for different epochs (ranged from 1 to 1000) is shown in Figure 5.4. As the number of epochs reaches 400, the error is not reduced any more as the model-fitting shows the over-fitting of the fine-tuning data and eventually the generalization capacity is reduced. Mixed-effect model: The performance of our proposed method is compared to the nonlinear mixed-effects model, which has been used extensively as a powerful growth hierarchical model 76 0.8 40 0.6 30 Training time (seconds) RMSE (cm) 0.4 20 0.2 10 0 0 0 200 400 600 800 1000 Number of epochs Figure 5.4: Effect of the number of epochs on the prediction error. The RMSE and the fine- tuning training time are plotted in solid blue and dashed red lines, respectively. The training time increases linearly with the number of epochs, while the RMSE rapidly decreases at the beginning but converges at around 400 epochs. over the decades [Brady et al., 2004, Eriksson et al., 2005, Sweeting and Thompson, 2012]. For the mixed-effects model, a basic form of the growth function is selected as: yi,j = α0 + (α1 + b1 )ti,j + (α2 + b2 )t2i,j + i,j , where yi,j and ti,j are the diameter and the associated time at the jth measurement of ith patient, b = [b1 , b2 ] is the random-effects terms and b ∼ N (0, Σb ), α = [α0 , · · · , α2 ] is the parameters vector, and i,j is the independent error term, i.e., i,j ∼ N (0, σw2 ). b and α are fitted to the data via the fminsearch function in MATLAB. 77 5.4 Test prediction results As it is shown in Table 5.2, a DBN with 300 nodes in both layers is selected for model testing. The absolute prediction error and relative prediction error for selected samples are shown in Table 5.3, and are compared with those from a mixed-effect model. A comparison of prediction results is shown in Figure 5.5. The proposed method outperforms the mix-effects method with a 65% reduction in RMSE. The overall prediction is relatively accurate because the relative error of each prediction only ranges from 2.3% to 4.3%, which is negligible. Table 5.3: The absolute and relative prediction errors of 6 testing samples under DBN and Mix-effects model. Sample ID RMSE (cm) Relative error* DBN Mix-effects DBN Mix-effects P1 0.224 0.393 4.3% 7.6% P2 0.181 0.645 3.1% 11.2% P3 0.168 0.535 2.8% 8.8% P4 0.147 0.655 2.3% 10.2% P5 0.197 0.406 2.9% 6.2% P6 0.165 0.494 3.1% 7.3% Mean value 0.180 0.521 3.1% 8.6% *The relative error is defined by the ratio between the absolute RMSE and the largest diameter in the objective IMDC, e.g, if the RMSE is 0.18 cm and the largest diameter is 5.8 cm, the relative error should be 0.18/5.8 = 3.1%. The deep learning model, which is implemented on the MATLAB, can be trained within 30 seconds on a PC with a 3.3 GHz 10-core CPU and a 64 GB RAM (Table 5.2). This period of time is short enough to provide insight in aiding clinicians to make surgical decision of AAAs. 78 P1 P2 P3 5 true true 4 true DL prediction 5 DL prediction DL prediction 4 ME prediction ME prediction ME prediction d (cm) d (cm) d (cm) 4 3 3 3 2 2 0 2 4 6 8 0 2 4 6 8 0 2 4 6 8 s (cm) s (cm) s (cm) P4 P5 P6 5 true true 4 DL prediction DL prediction 4 3.5 ME prediction ME prediction d (cm) d (cm) d (cm) 4 3 true 3 2.5 DL prediction 3 ME prediction 2 2 0 2 4 6 8 0 2 4 6 8 0 2 4 6 8 s (cm) s (cm) s (cm) Figure 5.5: The true value, the DBN model prediction and mixed-effects model prediction of IMDCs are shown in dotted dashed black (‘true’), solid red(‘DL prediction’), and dashed green lines(‘ME prediction’), respectively. s denotes the coordinate of location along the centerline of AAA and d represents the associated maximum diameter measured by inscribed sphere method. 5.5 Monte-Carlo cross-validation A Monte-Carlo cross-validation method is performed to show the robustness of the proposed deep learning model. As the first step, 13 sets of eligible testing data, i.e, {xi , yi }, are collected from 20 patients’ CT images. There are two criteria for choosing eligible data from a patient. The first criterion is that the eligible {xi , yi } should be the last set of data of the patient. The second criterion is that the number of raw CT images of the patient is at least four. As the second step, the cross-validation trials are independently performed 100 times under the deep structure of 2-layers DBN with 300 nodes on each layer. In each trial, we randomly choose 3 sets of test data out of the whole eligible dataset and leave the other sets of eligible data as training data to be used in the fine-tuning step. As a result of the Monte- Carlo cross-validation, 100 prediction errors (RMSEs) are independently collected, of which 79 the mean and the standard deviation are 0.196 cm and 0.051 cm. The standard deviation is so small relative to the mean that it guarantees the robustness of the proposed method. Moreover, as a comparison, the average testing RMSE (0.18 cm) falls into the range of the standard deviation of the cross-validation result, thus supporting the test results shown in Table.5.3. 5.6 Discussion This study utilized a physical G&R computational model combined with follow-up image data from 20 patients to predict the shape evolution of AAAs represented by IMDCs. To our knowledge, this is the first study that utilizes the deep learning technique to predict the shapes of AAAs in an evolutionary scheme based on a small dataset of follow-up images. The main difficulty in applying deep learning to predict AAA enlargement is the limited size of the training dataset, i.e, follow-up images of AAAs. In this study, we overcame this difficulty by proposing a work-flow, in which the DBN is pre-trained by massive in silico data. The accurate predictions demonstrate that deep learning holds promise for capturing evolutionary features of individual soft tissue with numerical simulations. It gives deep learning techniques a new application in biomedical engineering, other than surrogates [Liang et al., 2018], image segmentation [Lai, 2015], etc. Besides combining deep learning and AAA prediction, the proposed study also contributes by making fast and accurate predictions of AAA enlargement. Following Table 5.2, the model training time takes approximately 24 seconds, which is significantly faster than other statistical models [Do et al., 2019] [Zhang et al., 2019]. Additionally, the efficiency of data generation is significantly improved by the PCM approximation, which generates tens of thousands of IMDCs in 20 seconds. As a comparison, the G&R computational model takes 30 minutes to generate one set of data, meaning that it would take more than one year to 80 provide the same amount of training data. Additionally, due to the high complexity of physics in patient-specific predictions of AAA, even high-fidelity physical models cannot promise to make accurate predictions. To enable accurate predictions, in this study, we proposed a two-stage training approach: first, we pre-train our deep learning model with a computationally generated sizable dataset; second, we fine-tune the deep structure with the patient data for patient-specific predictions. In this study, the average prediction error is 0.180 cm, which is significantly small compared to the AAA diameters (3 cm to 8.5 cm). The results shown in Fig.5.5 and Table.5.3 also indicate that the proposed method indeed provides effective predictions, which outperform the clas- sical mixed-effect model [Brady et al., 2004, Eriksson et al., 2005, Sweeting and Thompson, 2012] by 65% in terms of the average relative error. Also, a Monte-Carlo cross-validation, which is a standard and widely recognized validating approach commonly used in machine learning studies, also provides similar prediction results compared with the test results, thus showing the effectiveness and robustness of the proposed study. Therefore, though the in sil- ico data are limited by the G&R which cannot capture all marginal situations, the two-stage DBN method still shows its promise in predicting patient-specific geometries. 81 Chapter 6 A PHYSICS-BASED MACHINE LEARNING APPROACH TO PREDICT ABDOMINAL AORTIC ANEURYSM ENLARGEMENT 6.1 Introduction Our motivation is to make an efficient patient-specific predictive tool, aiding in clinical treat- ment as a tool for predictive medicine. Nonetheless, a personalized simulation by its own is not capable of providing effective prediction. This is because the physical simulation is considered a deterministic approximation to the real patient data, which contains error compared with clinical measurements and cannot provide reliable confidence interval. Thus, given an accurate patient-specific simulation, developing a stochastic model, which com- bines both real-patient data and simulation outputs, is essential to enable patient-specific 82 prediction. Related methods that applied in vascular modeling include Bayesian calibration [Zhang et al., 2019], deep learning (Chapter 5) and machining learning [Liang et al., 2017]. Inspired by these studies, a novel physics-based machine learning approach is developed to extrapolate the limited set of clinical data (expensive data) and then further improve the accuracy of predictions by accommodating the patient-specific simulation outputs (cheap data). Other studies that combine the data from both experiments and simulations have shown the feasibility of prediction using the cokriging method [Fidkowski, 2014, Forrester, 2010, Kuya et al., 2011]. Therefore, cokriging method, as an efficient machine learning method, is perfectly qualified for this work because of it high efficiency and accuracy. In this section, we employ the cokriging approach, which is established in Section 3.3 for MFS, to enable another application, i.e., personalized prediction of AAA growth. The uncertainty analysis to the prediction is formulated by Equation 3.28. 6.2 Personalized prediction of AAA extension In this section, we employ cokriging as a data-driven predictive tool to predict AAA growth. In Section 4.2, we have developed an efficient parameter estimation algorithm that enables fast G&R simulations of individual patients. These simulations provide valuable information that supplements the clinical data at unmeasured spatial-temporal locations. For conve- nience, the G&R simulation results and clinical data of individual patients are denoted by ‘in silico data’ and ‘clinical data’, respectively. The in silico data and clinical data of each patient can be considered as the LF data and the HF data for training the cokriging based predictive tool. (1) (Ne ) For each patient, a set of Ne HF sample points with state variables Xe = {ξp , ..., ξp } is selected from the clinical observations, where the locations Xe are limited to the re- gion with available clinical observations. In addition, a set of Nc state variables Xc = 83 Physics-based Clinical Data G&R models cokriging Data-driven predictive tool of Prediction Multifidelity MFS-based arterial diseases surrogate parameter estimation algorithm In Silico Data Figure 6.1: Overall schematic flow of the proposed framework. A MFS is approximated from G&R models using a cokriging method. The in silico data of patient-specific simulations are obtained from both the MFS and the clinical data using a parameter estimation algorithm. The clinical and the in silico data are combined to predict the progression of arterial disease. (1) (N ) {ξˆs (θopt ), ..., ξˆs c (θopt )} can be sampled and substituted into the MFS of G&R to get the a set of Nc LF samples. In this case, the MFS is trained from the last LCB iteration during the parameter optimization, θopt represents the optimized patient-specific parameter, and ξˆs (θopt )) is the assimilated state variable estimated by Eq. (4.1). HF data and LF data are then combined to train a physics-based machine learning predictive tool of individual AAA expansion using cokriging (shown in Fig. 6.1). In the current setting of cokriging, the individualized simulation outputs can be regarded as the low-fidelity information that improves the limited high-fidelity data measured from real patients. According to the Eq. (3.24), if we let Ze (x) and Zc (x) represent the GPs of clinical data and the patient-specific in silico data, respectively, the GP of the correction term Zd (x) can be interpreted as an error term between them. All the hyperparameters can be trained using MLE. The final patient-specific prediction is expressed in the Eq. (3.27), and an explicit uncertainty analysis of the prediction error is expressed in Eq. (3.28). 84 55 55 50 50 45 45 d(mm) 40 d(mm) 40 35 35 30 30 25 25 15 20 25 30 35 40 15 20 25 30 35 40 s (mm) s (mm) (a) Patient 16 (b) Patient 10 Figure 6.2: Two examples that compare prediction results with the ground truth of AAA shape curve, where diameter (d) is plotted against the length along center line (s). The blue shades represent the 95% confidence interval of AAA growth prediction. The black circles and red crosses represent the ground truth of aneurysmal diameters that locate within and outside the confidence range, respectively. 6.3 Results of AAA expansion prediction Given best-fitted patient-specific parameters, we design a personalized data-driven predictive tool to predict the future growth of AAAs by combining both simulation outputs and clinical information (Section 6.2). For each patient whose number of CT images is greater than or equal to 3, we select the last CT image out of the follow-up CT images dataset as test data, and keep the remaining follow-up CT images as training data. Fig. 6.2 provides two examples of simulation results. Fig. 6.2(a) shows an example with a good prediction result, in which 9 out of 11 points are within 95% confidence interval. The average prediction error of diameter d = 0.24 mm is significantly smaller than 2σd = 0.96 mm, where σd represents the standard deviation associated with prediction. In Fig. 6.2(b), 7 out of 11 points locate within 95% confidence interval. Although the points on the right side of the AAA curve are not within the tolerance range, the overall prediction is acceptable. The average prediction error of diameter d = 1.28mm is less than 2σd = 1.78mm. 85 Table. 6.1 collects all prediction results and the associated uncertainty that can be directly compared with the ground truth. Note that Patient 2 only has 1 training CT scan that is not enough for parameter estimation, so we collect a total of 20 prediction results in the table. The computational time of each prediction is negligible (less than 1 minute) and has little effect on the overall computational time so the computational time of prediction is not recorded in the table. Table 6.1: Prediction results of personalized AAA expansion. Patient # Scans # Points Avg. Avg. σd Ground Predicted ID for within Prediction of truth of max. d a b c training tolerance error of d prediction max. d P01 2 9 0.29 0.45 51.26 51.10 P02 1 N/A N/A N/A N/A N/A P03 2 7 1.08 0.86 65.10 66.35 P04 4 11 0.47 0.67 51.02 50.35 P05 4 7 2.13 1.26 58.09 56.71 P06 4 9 0.83 1.22 34.74 37.33 P07 4 8 0.52 0.34 46.83 46.14 P08 4 2 3.23 0.78 55.09 50.76 P09 2 10 0.84 0.79 44.44 43.42 P10 3 7 1.28 0.89 46.91 48.49 P11 3 9 0.65 0.57 50.58 49.13 P12 2 8 0.09 0.17 46.61 46.35 P13 3 10 0.23 0.59 48.38 47.74 P14 2 0 5.05 0.51 51.54 46.48 P15 4 9 1.11 0.58 51.76 50.64 P16 2 9 0.24 0.48 48.07 49.21 P17 4 9 1.00 0.75 47.82 47.06 P18 3 9 0.12 0.66 51.50 52.70 P19 3 9 1.76 1.13 60.73 59.55 P20 3 6 1.84 0.86 51.11 48.66 P21 2 0 2.24 0.34 54.25 49.79 a The total number of data points is 11 for each patient. b Unit of d is mm. c σd represents the standard deviation of the prediction given by Eq. (3.28). We achieved satisfactory predictions for most patients. The mean value of average prediction error of arterial diameter among all patients is 1.39 mm, which is more accurate than the prediction in Chapter 5 where the average prediction error was 1.80 mm. The prediction 86 results of Patient 8, Patient 14 and Patient 21 are outside the tolerance range. It is partly because of the small number of scan images that limits the prediction performance. The longitudinal datasets of Patient 14 and Patient 21 only contain two CT images for MFS training, thereby reducing the accuracy of the parameter estimation and the prediction. Maximum diameter is an important indicator to assess the risk of AAAs [Paraskevas et al., 2010]. In this study, the average prediction error of the maximum diameter is 1.67 mm that is significantly small compared to AAA diameter (30 − 85 mm). In comparison, the prediction error of maximum diameter was 2.61 mm by using the patient-specific prediction from Akkoyun et al. [2020]; and Lee et al. [2018] can predict the maximum diameter within 2 mm error in 85% and 71% of patients at 12 and 24 months. In Chapter 5, we have developed a DBN with a two-stage training scheme that combines both a small amount of patient-specific data and a large amount of in silico data to predict AAA expansion. This DBN, however, only utilized the low-fidelity model, and the prediction capability of DBN can be improved by including a higher fidelity model. Accordingly, this study, which can provide more realistic high-fidelity G&R simulations, holds a great potential to be combined with our previous study of DBN to further improve the effectiveness of prediction. Therefore, we consider the proposed framework as a primary step towards a more sophisticated patient-specific simulation and prediction tool that will be capable of (1) integrating the information of physics-based modeling with big data, and (2) providing personalized clinical suggestions of arterial diseases within the clinically relevant time frame (1-2 hours) [Taylor et al., 1999]. 87 Chapter 7 A VIRTUAL PATIENT ENGINE FOR THE CREATION OF A MITRAL VALVE REGURGITATION PATIENT COHORT 7.1 Introduction Recent advances in biomechanics have led to an increasing trend of applying physics-based cardiovascular models to clinical-relevant applications, such as patient-specific simulations, testing medical treatments, and patient-specific disease progression prediction. However, limitations of computational models prevent the application of creating a cohort of in silico virtual patients as a testing environment for a specific clinical trial, i.e., In Silico Clinical Trial (ISCT). Main challenges in creating and analyzing a virtual patient population include • Computationally expensive 3D modeling 88 • Defining the salient characteristics (parameters) of each patient • Representing large patient variations • Obtaining anatomically plausible patients • Understanding the relationships between model parameters • Down selection to a specific targeted patient population • Etc. So as to overcome these challenges and generate in silico Virtual Patients Cohort (VPC) for a specific targeted patient population efficiently, we extend the proposed physics-based machine learning framework to create an innovative and general framework & strategy called Virtual Patient Engine (VPE). In particular, the VPE is built upon credential simulations of a parameterized physics-based model of the cardiovascular system. Next, the VPE uses machine learning (e.g., surrogate model) and statistical methods to efficiently down select cardiovascular simulations into a specific target population following empirical information of the targeted patients’ cohort such as geometrical ranges and mechanical characteristics. Finally, the VPE outputs a physics-based VPC containing a small set of physics-based Virtual Patients (VPs). Each VP includes a unique 3D finite element analysis of regurgitate mitral valve with desired characteristics. In this study, the VPE framework is demonstrated with an example problem of Secondary Mitral Regurgitation (SMR). This VPE employs a Mitral Valve Apparatus (MVA) model extracted from the Living Heart Human Model (LHHM) built on ABAQUS Explicit, and aim to generate VPs with physiologically feasible anatomy and a limited range of regurgitate MV opening size during systole. Here, the LHHM is a sophisticated platform supported by the Living Heart Project at Dassault Systèmes, which enables users to simulate different 89 heart diseases by adjusting features such as 3D geometries, coupled electrical-mechanical physics, hemodynamic parameters, and material parameters. This section is organized as followed. Section 7.2 and Section 7.3 describe the procedures for estimating the material parameter ranges and prescribing the parameterized geometry for mitral valves, respectively. An automatic simulation framework is developed in Section 7.4 to run simulations in batches. A Design of Experiment (DoE) is setup in Section 7.5 to generate MVA simulations. A machine learning powered VPE is developed in Section 7.6 to generate the VPC for SMR patients. Results and discussion are provided in Section 7.7 and 7.8, respectively. This chapter is the author’s research results during his internship at Dassault Systèmes. Ownership of the researches in this chapter fully resides with Dassault Systèmes. 7.2 Material parameter calibration Characterization of the values and ranges of mitral valve leaflet material parameters is a fundamental step towards the mitral valve simulations that covers large patient variations. To achieve this goal, we design the following workflow to calibrate the range of material parameters from experimental data: 1. Understand the stress-strain relationships (e.g. extensibility and stiffness) of leaflet tissue in both axial and circumferential directions from the uniaxial or biaxial tests in the literature [Grande-Allen et al., 2005, Pham and Sun, 2014]. 2. Convert the stress-strain relationships and their ranges into stress-strain data curves with uncertainty ranges. 3. Adjust the material parameters in the MVA constitutive model to match the stress- 90 strain data curve. 4. Collect the mean, lower limit and upper limit of material parameters. Two sets of material parameters are obtained by matching the data of two studies, including a biaxial testing from Pham and Sun [2014] and a uniaxial testing from Grande-Allen et al. [2005]. The summarized results are provided in Table 7.1. Table 7.1: Results of MVA material parameter estimation. Literature Attributes bf bs sf as Pham and Best fit 250 45 0.08 0.017 Sun [2014] range 250 (fixed) 45 (fixed) [0.04,0.15] [0.005,0.07] Grande-Allen Best fit 20 16 0.005 0.005 et al. [2005] range [6, 100] [4.8, 80] 0.005 (fixed) 0.005 (fixed) In this project, we perform VPE using the ranges of material parameters calibrated from the both references. The results show that the material parameters calibrated from Pham and Sun [2014] are too stiff to enable feasible mitral valve deformation. Thus, in this study, the range of material parameters is calibrated from data of Grande-Allen et al. [2005], which is marked in bold in Table 7.1. 7.3 Parameterization of mitral valve apparatus geom- etry In order to generate the virtual patient cohort that covers large inter-patient variability, each virtual patient in the cohort should contain a unique set of parameters and anatomy. According to Niederer et al. [2020], there are three strategies for generating virtual cohorts, including (1) 1:1 mapping virtual cohorts, (2) sampling from inferred distributions, and (3) random variation with acceptance criteria. The proposed VPE is based on the third strategy, i.e., random variation with acceptance criteria. This strategy requests to generate 91 VPs by randomly varying parameters and then include only these VPs that fall within the specific physiological bounds into the VPC. This means that the physics-based models used to generate VP simulations must cover various patient characteristics, such as the anatomic geometries and material characteristics, through changes in parameters. More details on VPC generation are discussed in later sections. In this section, we focus on the parameterization of the MVA model, which allows the geometrical parameters to capture anatomic geometries of mitral valve. A patient-specific MVA model is extracted from a patient-specific Living Heart Human Model (LHHM) as the basis for VPE creation. However, it is difficult to morph mitral valve geometry in a patient-specific MVA model to generate a large amount of VPs capturing inter- patient variety. Thus, in this section, we aim to extend the patient-specific MVA model into a more general parametric MVA model that can be easily morphed. A two-step strategy is implemented to achieve this goal. First, we analyze a large num- ber of mitral valve images and then summarize them into a couple of designs of mitral valve apparatus, including the annulus geometry, leaflet geometry, chordae insertion pattern, chor- dae origin pattern and papillary muscle locations. Second, these designs are automatically and efficiently translated into different geometries and meshes of MVA on the 3DExperience platform from Dassault Systèmes. This two-step parameterization is enabled by referring to anatomical literature, especially Oliveira et al. [2020]. Due to confidentiality issues, the detailed design is not provided in this dissertation. A demo simulation is performed and illustrated in Figure 7.1. 92 Figure 7.1: A demo simulation of parameterized mitral valve apparatus model. . 7.4 Automatic simulation framework An automatic simulation framework is developed to efficiently generate a large number of simulations to create a surrogate model and a virtual patient cohort. This framework can be executed on a server installed with ABAQUS or on the 3DS Experience platform. The process of the automatic simulation framework is described as follows: 1. Generate geometries and mesh files for parameterized MVA models with different input parameters. 2. Assemble mesh files and other parameters to produce simulation input files of ABAQUS. 3. Submit simulation input files to a server installed with ABAQUS to run simulations in batches and obtain simulation output files. 4. Post-process the simulation output files to obtain output target quantities. In this simulation framework, Python scripts automatically assemble meshes, produce sim- ulation input files, submit input files to server, and perform post-processing. Thus, this 93 Figure 7.2: A diagram of the automatic simulation framework for parameterized MVA model. . working is considered to be fully automatic, which provides a basis to create the VPE. A di- agram of the automatic simulation framework is shown in Figure 7.2. The input and output parameters will be automatically saved in a table after completing executing the simulation framework. As indicated in Figure 7.2, post-processing scripts are developed to extract sev- eral target quantities from the simulation output files. In this study, four main anatomic quantities are selected as target quantities, including anatomic mitral valve regurgitate ori- fice area, tenting angle of leaflets, tenting length of leaflets and coaptation length of leaflets. These geometrical quantities are selected because they are salient features for surgeons to capture the individual characteristics and severity of the mitral valve orifice area. 7.5 Design of experiment setup We created an automatic simulation framework for parameterized MVA model. Next, to cover the inter-patient variety with various physiology, anatomy, mechanics, etc., we need to produce a large number of simulations of VPs with different values of input parameters. Specifically, we need to carefully select input parameters to satisfy four requirements. First, most of the parameters that strongly affect target outputs are included. Second, the def- inition of input parameters should be specified by leveraging available experimental data, clinical information, and modeling experiences to ensure that the simulations are physiolog- ically feasible. Third, the ranges of input parameters should be determined by the literature related to the targeted patient cohort as well as engineering experiences. Last, the values of 94 less important parameters are fixed with population mean of the targeted population after sensitivity analysis. In this study, we select 10 parameters to propagate inter-patient variations into VPs, includ- ing 8 geometrical parameters and 2 material parameters, as provided in Table 7.2. Their ranges are determined by anatomic literature of secondary mitral valve regurgitation [Kun- zelman et al., 1994, Oliveira et al., 2020, Veronesi et al., 2008]. Note that the LLR, leaflet length ration, is define by Lengthof Leaf let(u) BaslineAP diameter LLR = × . (7.1) BaselineLengthof Leaf let(u) AP diameter Table 7.2: List of input parameters and their ranges. List of Range Definition of parameters parameters APD [32.4, 45.2] mm AP diameter of annulus ADR [0.94,1.24] Ratio between AP and ALPM diameter of annulus LLR [0.8,1.2] Leaflet length ratio AH [4,7.2] mm Annulus height APM [27.5,44.5] mm The vertical distance between annulus plane and anterior papillary muscle PPM [28.7,46.5] mm The vertical distance between annulus plane and posterior papillary muscle Phi [85, 135] degree The angular position of papillary muscles relative to mitral valve RP [2,8] mm The radius of C-Shape chordae origin points Pre-stretch [-0.1,0.1] Chordae pre-stretch Bf [6, 100] Leaflet tissue stiffness In order to reduce the computational time of creating a surrogate model and to generate a VPC within a small amount of time, e.g., one day, the input parameters of MVA simulations should be selected from the 10-dimensional parameter space through an optimized strategy. In this study, we employed the Latin Hypercube Sampling (LHS) method which has the advantage in sampling from high-dimensional distributions. Typically, the number of samples required to train the Kriging surrogate model (Section 3.3.1) should be equal to 10∼20 times 95 the number of dimensions, thus we use LHS to select 250 samples, of which 200 samples are used for and 50 samples are used for testing. Next, MVA simulations are performed on 14-core computational nodes on an LSF serve. Each simulation takes between 25 to 40 minutes. Thus, if all 250 simulations were performed in a pipeline, it would take about 140 hours to complete all simulations. Fortunately, we submitted simulations to the LSF server altogether and multiple simulations ran in parallel. Therefore, in practice, it takes 10 hours to complete all simulations. Moreover, the results show the robustness of the proposed the simulation framework. All simulations and the post-processing are fully automated with a 100% successful rate. A demonstration of the simulation results is shown in Table 7.3, where 12 out of 250 samples are provided in the rows, and the input parameters and the target output (orifice area) are shown in the columns. Table 7.3: A demonstration of LHS samples with the values of input parameters and targeted outputs. Input parameters Output In- APD ADR LLR AH Pre- APM PPM Phi RP Bf Orifice dex stretch area (mm2 ) 1 40.9 0.976 1.04 6.92 -0.0211 34.98 32.24 95.8 6.05 53.39 143.8 2 41.0 0.987 1.09 4.58 0.0283 31.27 45.45 125.2 3.66 55.05 54.3 3 39.1 1.105 1.16 6.12 -0.0217 38.08 40.49 108.2 2.67 8.81 1.3 4 32.9 1.10 1.19 5.43 0.0888 33.01 41.00 88.9 7.77 34.01 10.9 5 42.3 1.167 1.19 6.63 0.0414 35.76 43.33 89.7 5.44 68.31 120.2 6 42.6 0.970 0.99 4.03 0.0925 34.42 35.90 123.5 3.17 9.36 24.8 7 43.1 1.163 0.994 5.36 0.0871 33.17 33.09 90.7 3.37 61.20 191.2 8 37.9 1.087 0.81 5.25 -0.0651 27.64 35.27 134.1 2.38 50.63 231.5 9 37.6 1.071 0.91 4.13 0.0316 33.99 46.13 109.6 6.13 12.43 66.5 10 37.1 1.108 1.17 6.90 -0.0839 40.03 32.94 120.8 7.69 37.15 2.4 11 41.6 1.236 1.07 7.06 0.0981 37.02 36.16 102.3 2.78 21.05 46.3 12 38.0 1.039 0.91 4.66 0.0518 43.14 29.24 89.4 5.96 10.87 95.8 ... ... ... ... ... ... ... ... ... ... ... ... 96 Figure 7.3: Simulation outputs vs. surrogate outputs of different surrogate models. . 7.6 Virtual patient engine In this section, we describe the workflow for creating a VPE from DoE simulations. 7.6.1 Create surrogate model In order to create a surrogate model, 250 samples from DoE are divided into two sub- datasets, including a training dataset and a testing dataset. In the training of the surrogate, we further split the training dataset into different folds of data for k times following the rules of k-fold cross-validation, and then implement a Bayesian auto-tuning algorithm to find the best-estimated hyper-parameters of surrogate. Next, we perform model testing to the surrogates by matching the simulation results from the original testing dataset with the approximated simulation results from the surrogate model (surrogate outputs). We train and test different surrogates created by different types of machine learning methods, such as kriging, kernel ridge and XGBoost, from a Python package, scikit-learn [Buitinck et al., 2013]. The model testing results are listed in Table 7.4. The surrogate with the best performance, i.e. the kriging with Matern kernel and hyper- parameters auto-tuned by Bayesian optimization, gives a mean absolute error of 12.9 mm2 of orifice area and an R2 of 0.959. 97 Table 7.4: Performance of different surrogate approaches in predicting anatomic orifice area (mm2 ). Method Kernel Mean absolute Mean squared R2 error (mm2 ) error Kriging Default 13.4 421 0.946 Kriging Matern (1000,1) 13.1 329 0.958 Kriging Rational Quadratic (10,0.005) 14.7 406 0.950 Kriging Matern (auto-tuned) 12.9 322 0.959 Kernel ridge Kernel Ridge (0.001) 27.9 1388 0.821 XGBoost Default 23.8 993 0.862 XGBoost Auto-tuned 23.6 949 0.868 Figure 7.4: Simulation outputs vs. surrogate output of other target quantities, including tenting angle of the anterior leaflet, the tenting angle of posterior leaflet and the tenting height. . Figure 7.3 provides four plots comparing the performance of four surrogate approaches in predicting orifice area, including an XGBoost, a kernel ridge, a kriging with the default setup and a kriging with the hyper-parameters auto-tuned. In addition, Figure 7.4 shows modeling testing results to the other target quantities, including the tenting angle of anterior leaflet, the tenting angle of posterior leaflet and the tenting height. Note that all surrogate models shown in Figure 7.4 are trained by kriging. 7.6.2 Generate surrogate-based VPC Given a well-trained surrogate, we are able to generate a large amount of surrogate-based VPs with various characteristics. In order to constrain the variation of simulation results to 98 Figure 7.5: Histogram of orifice area for surrogate-based VPs. . a reasonable range and produce VPs with feasible physiological and mechanical properties, the input parameters are assumed to be independent and following the Gaussian distribution with a mean of 0.5 and a standard deviation of 0.33 after normalizing each parameter into the range [0, 1]. Next, Monte-Carlo sampling is implemented to sample 100,000 sets of input parameters as input parameters of the surrogate-based VPC. Next, those input parameters are substituted into the surrogate to generate the surrogate outputs, i.e., 100,000 outputs of orifice area simulated by the surrogate. Figure 7.5 gives a histogram of orifice area for all surrogate-based VPs, showing that most VPs have orifice areas between 0 ∼ 150 mm2 which is consistent with the reference values of mitral valve regurgitation. 7.6.3 Down select VPC for a targeted population As indicated by literature [Baumgartner et al., 2017, Otto et al., 2021], patients with re- gurgitate orifice areas between 30 and 100 mm2 may exhibit moderate to severe mitral 99 Figure 7.6: Histogram of orifice area for surrogate-based VPs and the down selection crite- rion. . regurgitation and can be considered for repair through clipping. Note that this criterion is a simple and temporary inclusion criterion designed for demonstrating the VPE. It does not guarantee a VPC that the physiology and mechanics features are exactly same as a real pa- tient cohort. In the future, the inclusion criteria can be further refined as additional evidence becomes available. Figure 7.6 illustrates how inclusion criterion performs down selection to the surrogate-based samples. As indicated in Figure 7.7, about 40,000 out of 100,000 Monte-Carlo samples are selected in the down selection. Although the input parameters associated with the Monte-Carlo samples are specified to be independent of each other, correlations may arise in the samples after performing the down selection. For example, Figure 7.7 gives a demo with 200 Monte- Carlo samples, in which the included samples marked by red show a clear trend of diagonal distribution. More specifically, the samples in the upper-left region are likely to correspond to the MVs with extremely large orifices, while the samples on the lower-right region are likely to correspond to the MVs with closed orifices. 100 Figure 7.7: A demonstration figure with 200 Monte-Carlo samples displayed on a two- dimensional parameter space of Bf vs. leaflet ratio (LLR). Included samples are indicated by red dots, and the excluded samples are indicated by blue circles. . 7.6.4 Find representative VPs using a clustering algorithm As aforementioned, the primary goal of VPE is to select some sets of input parameters, which could cover the inter-patient variety of desired characteristics after being imposed into the physics-based model. This requires us to select a small number of samples (50 samples in this study) that are representative of different variations of the target population. In contrast, there are around 40,000 surrogate-based samples included through the down selection. Thus, another down sampling should be performed to further select 50 representative VP samples from 40,000 down selected samples, and we can run physics-based simulations with the input parameter sets of these 50 VP samples. To achieve this goal, in this study, a computationally efficient machine learning method, i.e. clustering, is used. Clustering is an unsupervised machine learning task which automatically 101 classify data points according to their natural patterns. A large number of clustering meth- ods, such as k-mean, mean shift and spectral clustering, can be easily to be implemented using a Python package, scikit-learn [Buitinck et al., 2013]. However, most of these clustering methods cannot provide good results due to the ‘high-dimensional curse’. Fortunately, after studying many alternative approaches, we find that the affinity propaga- tion approach is able to address the ‘high-dimensional curse’. It can select cluster centers that evenly cover the whole parameter space. Along with another method, a factor-2 approx- imation method, which allows to extract 50 VP samples from cluster centers by maximizing relative distance between cluster centers, the affinity propagation approach produces input parameters of 50 VP samples that cover a large inter-patient variety and many worst-case scenarios. Figure 7.8 shows a 3D representation of parameter space, in which 50 VP samples are represented by the red dots and a sub-group of the included Monte-Carlo samples are represented by green circles. This figure is shown in a 3D parameter space (APD vs. LLR vs. Bf) projected from a 10-dimensional parameter space. 7.7 Results: physics-based VPC and statistical analy- sis After collecting all the input parameter values of 50 VP samples, we could substitute these values into the physics-based parameterized MVA model to obtain the physics-based VPC with desired characteristics. Some results for physics-based VPC are shown in Figure 7.9, in which we can observe top- down and A2-P2 cut views of mitral valve. Figure 7.9 shows the main target quantities including the coaptation length, tenting angle and orifice area. In general, we see a long valley-shaped gap or double orifice. Figure 7.9 also indicates the values of gap size and 102 Figure 7.8: Included Monte-Carlo samples and 50 VP samples in a 3D parameter space. . 103 Figure 7.9: Top-down view and A2P2 cut view of 12 VPs in the VPC. orifice area for VP 9 which is about the maximum. As a first sanity check on the plausibility of this cohort, we see that our outputs are fairly consistent with those reported in the literature for patients with SMR. For more details, Table 7 provides the target quantities extracted from the physics-based VPC. Figure 7.10 provides a histogram of the orifice area extracted from the physics-based VPC. As a comparison, Nishimura et al. [2014] and Baumgartner et al. [2017] cite that an EROA ≥ 20mm2 is an indicator of severe SMR, and Oliveira et al. [2020] reports an average tenting height of 9.7 ± 3.2 mm and average posterior leaflet tenting angles of 44.4 ± 11.9◦ for severe MR patients. 7.8 Discussion In summary, by integrating FEA and Machine learning, this study creates a novel general framework and strategy of VPE which holds the potential to efficiently reproduce the VPC 104 Figure 7.10: A histogram of the orifice area extracted from 50 virtual patients in the physics- based VPC. of the targeted population with a variety of characteristics. This study is demonstrated by generating the VPs with SMR. It takes 10 hours to run 250 finite element simulations of MVA to provide input of VPE, and takes around 3 hours to execute VPE to generate 50 VP cohorts of the targeted population. The proposed VPE holds promise for many applications. Firstly and most importantly, the VPC generated by VPE can be used to provide a testing environment for medical treatments as an ISCT. Secondly, the VPE allows to prescribe the inter-patient variety among cohort members, which enables the analysis between the parameters, disease characteristics and the performance of clinical trials. Thirdly, the VPE enables a data augmentation workflow that efficiently reproduces a massive amount (N > 100,000) of surrogate-based simulation outputs. These surrogate-based simulations allow to analyze the correlations among dif- ferent factors of diseases such as material properties of biological tissue, anatomy, clinical measurable information. Fourthly, the VPE can efficiently correlate input parameters and the simulation outputs so that it holds the potential to propagate uncertainties of inputs, clinical measurements and modeling error into an estimation of targeted simulation outputs. Fifthly, VPE can be further extended to generate big data to train deep learning predictive models to better capture the disease progression. Lastly, the surrogate model in VPE con- 105 tributes by providing a fast inference between input parameters (e.g., material properties and geometry of organ) and measurable and unmeasurable clinical outputs. Though a lot of progress has been made by the proposed VPC, we are keenly aware that limitations exist in this preliminary study. First, it is not sufficient to claim it’s a plausible cohort. And as such one of our next tasks will be focusing on how to do a thorough assessment of the VPC in order to confidently state that it is representative of the real world. Second, even if the VPC is representative, it’s quite possible that the VPC corresponds to only a subset of the overall population. In the future, more input parameters may need to be defined and more variation in the existing parameters may be necessary to broaden the scope of the cohort to represent a larger sample of the overall MVR population. We will use the same VPE framework to adapt to the DoE of physics-based simulations with more varied parameters. Third, this study focuses on a simplified MVA model which is too simple to represent the human physiology of mitral valve with complex interactions with blood circulation, left ventricle, etc. In the future, models with higher fidelity should be introduced to further demonstrate that the VPC can be applied for the clinical practice of ISCTs. Last, in this study, we lack sufficient patient data to support the parameterization of mitral valve geometry, ranges of parameters, and correlation among parameters. In the future, if rich clinical data are collected, we could create multiple one-to-one mapping patient-specific models to build better-parameterized models and use patient-specific data to fine-tune the values, ranges, and correlations among parameters by using statistical methods such as Bayesian calibration. All of these works remain to be done before using the cohort in a predictive iSCT. 106 BIBLIOGRAPHY 107 BIBLIOGRAPHY F. Abraham, M. Behr, and M. Heinkenschloss. Shape optimization in unsteady blood flow: A numerical study of non-newtonian effects. Computer methods in biomechanics and biomedical engineering, 8:201–12, 07 2005. E. Akkoyun, S. Kwon, A. Acar, W. Lee, and S. Baek. Predicting abdominal aortic aneurysm growth using patient-oriented growth models with two-step bayesian inference. Computers in Biology and Medicine, 117:103620, 02 2020. M. S. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, and G. N. Wells. The fenics project version 1.5. Archive of Numerical Software, 3(100), 2015. D. Ambrosi, M. B. Amar, C. J. Cyron, A. DeSimone, A. Goriely, J. D. Humphrey, and E. Kuhl. Growth and remodelling of living tissues: perspectives, challenges and opportu- nities. Journal of The Royal Society Interface, 16:20190233, 2019. S. Baek, K. R. Rajagopal, and J. D. Humphrey. Competition between radial expansion and thickening in the enlargement of an intracranial saccular aneurysm. Journal of Elasticity, 80:13–31, 2005. S. Baek, K. R. Rajagopal, and J. D. Humphrey. A theoretical model of enlarging intracranial fusiform aneurysms. Journal of Biomechanical Engineering, 128(1):142, 2006. S. Baek, A. Valentin, and J. Humphrey. Biochemomechanics of cerebral vasospasm and its resolution: Ii. constitutive relations and model simulations. Annals of biomedical engi- neering, 35:1498–509, 10 2007. S. Balocco, O. Camara, E. Vivas, S. Teresa, L. Guimaraens, H. Gratama van Andel, C. Ma- joie, J. Pozo, B. Bijnens, and A. Frangi. Feasibility of estimating regional mechanical properties of cerebral aneurysms in vivo. Medical physics, 37:1689–706, 04 2010. H. Baumgartner, V. Falk, J. J. Bax, M. De Bonis, C. Hamm, P. J. Holm, B. Iung, P. Lan- cellotti, E. Lansac, D. Rodriguez Muñoz, R. Rosenhek, J. Sjögren, P. Tornos Mas, A. Va- hanian, T. Walther, O. Wendler, S. Windecker, J. L. Zamorano, and ESC Scientific Docu- ment Group. 2017 ESC/EACTS Guidelines for the management of valvular heart disease. European Heart Journal, 38(36):2739–2791, 08 2017. Y. Bengio. Learning deep architectures for AI. Foundations and trends® in Machine Learn- ing, 2(1):1–127, 2009. Y. Bengio, P. Simard, and P. Frasconi. Learning long-term dependencies with gradient descent is difficult. IEEE transactions on neural networks, 5(2):157–166, 1994. 108 R. H. Bhak, M. Wininger, G. R. Johnson, F. A. Lederle, L. M. Messina, D. J. Ballard, S. E. Wilson, and ADMS Group. Factors associated with small abdominal aortic aneurysm expansion rate. JAMA Surgery, 150(1):44–50, 01 2015. A. Booker, J. Dennis, P. Frank, D. Serafini, V. Torczon, and M. Trosset. A rigorous frame- work for optimization of expensive functions by surrogates. Structural Optimization, 17, 09 1998. A. Borghi, N. Wood, and X. Xu. Fluid-solid interaction simulation of flow and stress pat- tern in thoracoabdominal aneurysms: A patient specific study. Journal of Fluids and Structures, 24:270–280, 02 2008. A. R. Brady, S. G. Thompson, F. G. R. Fowkes, R. M. Greenhalgh, J. T. Powell, and UK Small Aneurysm Trial Participants. Abdominal aortic aneurysm expansion risk factors and time intervals for surveillance. Circulation, 110(1):16–21, 2004. P. Brown, D. Zelt, and B. Sobolev. The risk of rupture in untreated aneurysms: The impact of size, gender, and expansion rate. Journal of vascular surgery : official publication, the Society for Vascular Surgery and International Society for Cardiovascular Surgery, North American Chapter, 37:280–4, 03 2003. L. Buitinck, G. Louppe, M. Blondel, F. Pedregosa, A. Mueller, O. Grisel, V. Niculae, P. Pret- tenhofer, A. Gramfort, J. Grobler, R. Layton, J. VanderPlas, A. Joly, B. Holt, and G. Varo- quaux. API design for machine learning software: experiences from the scikit-learn project. In ECML PKDD Workshop: Languages for Data Mining and Machine Learning, pages 108–122, 2013. W. A. Cappeller, H. Engelmann, S. Blechschmidt, M. Wild, and L. Lauterjung. Possible objectification of a critical maximum diameter for elective surgery in abdominal aortic aneurysms based on one- and three-dimensional ratios. The Journal of cardiovascular surgery, 38:623–8, 01 1998. A. Caroli, S. Manini, L. Antiga, K. Passera, B. Ene-Iordache, S. Rota, G. Remuzzi, A. Bode, J. Leermakers, F. van de Vosse, R. Vanholder, M. Malovrh, J. Tordoir, and A. Remuzzi. Validation of a patient-specific hemodynamic computational model for surgical planning of vascular access in hemodialysis patients. Kidney international, 84, 05 2013. Y. Cheng, F. Wang, P. Zhang, and J. Hu. Risk Prediction with Electronic Health Records: A Deep Learning Approach. In Proceedings of the 2016 SIAM International Conference on Data Mining, 2016. J. Choi, J. Lee, and S. Oh. Swarm intelligence for achieving the global maximum using spatio-temporal gaussian processes. In 2008 American Control Conference, pages 135– 140, 06 2008. E. Choke, G. Cockerill, W. Wilson, S. Sayed, J. Dawson, I. M. Loftus, and M. Thompson. A review of biological factors implicated in abdominal aortic aneurysm rupture. European journal of vascular and endovascular surgery : the official journal of the European Society for Vascular Surgery, 30:227–44, 10 2005. 109 R. Collobert and J. Weston. A unified architecture for natural language processing: Deep neural networks with multitask learning. In Proceedings of the 25th international confer- ence on Machine learning, pages 160–167. ACM, 2008. M. Colombo, M. Bologna, M. Garbey, S. Berceli, Y. He, J. Rodriguez, F. Migliavacca, and C. Chiastra. Computing patient-specific hemodynamics in stented femoral artery models obtained from computed tomography using a validated 3d reconstruction method. Medical engineering & physics, 10 2019. I. Couckuyt, T. Dhaene, and P. Demeester. ooDACE toolbox: A flexible object-oriented Kriging implementation. Journal of Machine Learning Research, 15:3183–3186, 2014. C. Cyron, R. Aydin, and J. D. Humphrey. A homogenized constrained mixture (and mechan- ical analog) model for growth and remodeling of soft tissue. Biomechanics and Modeling in Mechanobiology, 15:1389–1403, 2016. H. N. Do, A. Ijaz, H. Gharahi, B. Zambrano, J. Choi, W. Lee, and S. Baek. Prediction of abdominal aortic aneurysm growth using dynamical gaussian process implicit surface. IEEE Transactions on Biomedical Engineering, 66(3):609–622, 2019. V. Dubourg, B. Sudret, and J. Bourinet. Reliability-based design optimization using kriging surrogates and subset simulation. Structural and Multidisciplinary Optimization, 44:673– 690, 11 2011. P. Eriksson, S. Jormsjö-Pettersson, A. R. Brady, H. Deguchi, A. Hamsten, and J. T. Powell. Genotype–phenotype relationships in an investigation of the role of proteases in abdominal aortic aneurysm expansion. British Journal of Surgery, 92(11):1372–1376, 2005. N. Famaey, J. Vastmans, H. Fehervary, L. Maes, E. Vanderveken, F. Rega, S. J. Mousavi, and S. Avril. Numerical simulation of arterial remodeling in pulmonary autografts. Journal of Applied Mathematics and Mechanics, pages 2239–257, 2018. M. Farsad, S. Zeinali-Davarani, J. Choi, and S. Baek. Computational growth and remod- eling of abdominal aortic aneurysms constrained by the spine. Journal of Biomechanical Engineering, 137(9):091008, 2015. G. Fernandez, C. Park, N. Kim, and R. Haftka. Review of multi-fidelity models. arXiv:1609.07196, 03 2017. K. J. Fidkowski. Quantifying uncertainties in radiation hydrodynamics models. 2014. A. Forrester. Black-box calibration for complex-system simulation. Philosophical transac- tions. Series A, Mathematical, physical, and engineering sciences, 368:3567–79, 08 2010. A. Forrester, A. Sobester, and A. Keane. Multi-fidelity optimization via surrogate modelling. Proc. R. Soc. A, 463:3251–3269, 12 2007. K. Fukushima. Neocognitron: A self-organizing neural network model for a mechanism of pattern recognition unaffected by shift in position. Biological cybernetics, 36(4):193–202, 1980. 110 H. Gharahi, B. A. Zambrano, C. Lim, J. Choi, W. Lee, and S. Baek. On growth measurements of abdominal aortic aneurysms using maximally inscribed spheres. Medical Engineering & Physics, 37(7):683–691, 2015. A. Ghorpade and B. Baxter. Biochemistry and molecular regulation of matrix macro- molecules in abdominal aortic aneurysms. Annals of the New York Academy of Sciences, 800:138–50, 12 1996. X. Glorot, A. Bordes, and Y. Bengio. Deep sparse rectifier neural networks. In Proceedings of the Fourteenth International Conference on Artificial Intelligence and Statistics, pages 315–323, 2011. K. Grande-Allen, J. Barber, K. Klatka, P. Houghtaling, I. Vesely, C. Moravec, and P. Mc- Carthy. Mitral valve stiffening in end-stage heart failure: evidence of an organic con- tribution to functional mitral regurgitation. The Journal of thoracic and cardiovascular surgery, 130:783–90, 10 2005. C. He and M. Roach. The composition and mechanical properties of abdominal aortic aneurysms. Journal of vascular surgery : official publication, the Society for Vascular Surgery and International Society for Cardiovascular Surgery, North American Chapter, 20:6–13, 08 1994. X. He, R. S. Zemel, and M. Carreira-Perpiñán. Multiscale conditional random fields for image labeling. In Computer vision and pattern recognition, 2004. CVPR 2004. Proceedings of the 2004 IEEE computer society conference on, volume 2, pages II–695. IEEE, 2004. G. E. Hinton. Training products of experts by minimizing contrastive divergence. Neural computation, 14(8):1771–1800, 2002. G. E. Hinton. Learning multiple layers of representation. Trends in cognitive sciences, 11 (10):428–434, 2007. G. E. Hinton. A practical guide to training restricted Boltzmann machines. Momentum, 9 (1):926, 2010. G. E. Hinton, T. J. Sejnowski, and D. H. Ackley. Boltzmann machines: Constraint satisfac- tion networks that learn. Carnegie-Mellon University, Department of Computer Science Pittsburgh, PA, 1984. G. E. Hinton, S. Osindero, and Y. Teh. A fast learning algorithm for deep belief nets. Neural computation, 18(7):1527–1554, 2006. G. E. Hinton, N. Srivastava, A. Krizhevsky, I. Sutskever, and R. R. Salakhutdinov. Improving neural networks by preventing co-adaptation of feature detectors. CoRR, abs/1207.0580, 2012. K. Hirata, T. Nakaura, M. Nakagawa, M. Kidoh, S. Oda, D. Utsunomiya, and Y. Yamashita. Machine learning to predict the rapid growth of small abdominal aortic aneurysm. Journal of Computer Assisted Tomography, 44:37–42, 2020. 111 G. Holzapfel, T. Gasser, and R. Ogden. A new constitutive framework for arterial wall mechanics and a comparative study of material models. Journal of Elasticity, 61:1–48, 04 2012. W. Huang, G. Song, H. Hong, and K. Xie. Deep architecture for traffic flow prediction: deep belief networks with multitask learning. IEEE Transactions on Intelligent Transportation Systems, 15(5):2191–2201, 2014. J. D. Humphrey and K. R. Rajagopal. A constrained mixture model for growth and re- modeling of soft tissues. Mathematical Models and Methods in Applied Sciences, 12(3): 407–430, 2002. L. M. Itu, P. Sharma, T. Passerini, A. Kamen, C. Suciu, and D. Comaniciu. A param- eter estimation framework for patient-specific hemodynamic computations. Journal of Computational Physics, 281, 01 2015. Z. Jiang, H. N. Do, J. Choi, W. Lee, and S. Baek. A deep learning approach to predict abdominal aortic aneurysm expansion using longitudinal data. Frontiers in Physics, 7: 235, 2020. D. R. Jones, M. Schonlau, and W. J. Welch. Efficient global optimization of expensive black-box functions. Journal of Global Optimization, 13(4):455–492, Dec 1998. M. Jordanski, M. Radovic, Z. Milosevic, N. Filipovic, and Z. Obradovic. Machine learning approach for predicting wall shear distribution for abdominal aortic aneurysm and carotid bifurcation models. IEEE Journal of Biomedical and Health Informatics, 22(2):537–544, 2018. A. Kahn, J. Marzat, H. Piet-Lahanier, and M. Kieffer. Global extremum seeking by kriging with a multi-agent system. IFAC-PapersOnLine, 48(28):526 – 531, 2015. 17th IFAC Symposium on System Identification SYSID 2015. S. Karnik, B. Brooke, A. Bayes-Genis, L. Sorensen, J. Wythe, R. Schwartz, M. Keating, and D. Li. A critical role for elastin signaling in vascular morphogenesis and disease. Development (Cambridge, England), 130:411–23, 02 2003. M. Kennedy and A. O’Hagan. Predicting the output from a complex computer code when fast approximations are available. Biometrika, 87, 10 1998. D. Klepach, L. C. Lee, J. F. Wenk, M. B. Ratcliffe, T. I. Zohdi, J. A. Navia, G. S. Kassab, E. Kuhl, and J. M. Guccione. Growth and remodeling of the left ventricle: A case study of myocardial infarction and surgical ventricular restoration. Mechanics Research Com- munications, 42:134–141, 2012. A. Klink, F. Hyafil, J. Rudd, P. Faries, V. Fuster, Z. Mallat, O. Meilhac, W. Mulder, J. Michel, F. Ramirez, G. Storm, R. Thompson, I. Turnbull, J. Egido, J. Martı́n-Ventura, C. Zaragoza, D. Letourneur, and Z. Fayad. Diagnostic and therapeutic strategies for small abdominal aortic aneurysms. Nature Reviews Cardiology, 8(6):338–347, 2011. 112 E. Kuhl, R. Maas, G. Himpel, and A. Menzel. Computational modeling of arterial wall growth : Aaattempts towards patient-specific simulations based on computer tomography. Biomechanics and modeling in mechanobiology, 6:321–31, 10 2007. K. Kunzelman, R. P. Cochran, E. Verrier, and R. Eberhart. Anatomic basis for mitral valve modeling. The Journal of heart valve disease, 3:491–6, 10 1994. Y. Kuya, K. Takeda, X. Zhang, and A. Forrester. Multifidelity surrogate modeling of exper- imental and computational aerodynamic data sets. Aiaa Journal - AIAA J, 49:289–298, 02 2011. S. T. Kwon, W. Burek, A. C. Dupay, M. Farsad, E. A. Park, W. Lee, and S. Baek. Interaction of expanding abdominal aortic aneurysm with surrounding tissue: Retrospective CT image studies. Journal of Nature and Science, 1:e150, 2015. S. K. Kyriacou, C. Schwab, and J. D. Humphrey. Finite element analysis of nonlinear orthotropic hyperelastic membranes. Computational Mechanics, 18:269–278, 08 1996. M. Lai. Deep learning for medical image segmentation. arXiv:1505.02000, 2015. B. L. Langille. Arterial remodeling: Relation to hemodynamics. Canadian journal of physi- ology and pharmacology, 74:834–41, 08 1996. J. D. Laubrie, J. S. Mousavi, and S. Avril. A new finite element shell model for arterial growth and remodeling after stent implantation. International Journal for Numerical Methods in Biomedical Engineering, 36, 2019. Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner. Gradient-based learning applied to docu- ment recognition. Proceedings of the IEEE, 86(11):2278–2324, 1998a. Y. LeCun, C. Cortes, and C. J. C. Burges. The MNIST database of handwritten digits, 1998b. J. Lee, Y. Bahri, R. Novak, S. Schoenholz, J. Pennington, and J. Sohl-Dickstein. Deep Neural Networks as Gaussian Processes. arXiv e-prints, art. arXiv:1711.00165, 10 2017a. R. Lee, A. Jones, I. Cassimjee, and A. Handa. International opinion on priorities in research for small abdominal aortic aneurysms and the potential path for research to impact clinical management. International Journal of Cardiology, 245:253–255, 10 2017b. R. Lee, D. Jarchi, R. Perera, A. Jones, I. Cassimjee, A. Handa, D. Clifton, K. Bellamkonda, F. Woodgate, N. Killough, N. Maistry, A. Chandrashekar, C. R. Darby, A. Halliday, L. J. Hands, P. Lintott, T. R. Magee, A. Northeast, J. Perkins, and E. Sideso. Applied machine learning for the prediction of growth of abdominal aortic aneurysm in humans. EJVES Short Reports, 39, 05 2018. D. Li, B. Brooke, E. Davis, R. Mecham, L. Sorensen, B. Boak, E. Eichwald, and M. Keating. Elastin is an essential determinant of arterial morphogenesis. Nature, 393:276–80, 06 1998. 113 T. Li and X. L. Yang. An efficient uniform design for kriging-based response surface method and its application. Computers and Geotechnics, 109:12–22, 01 2019. L. Liang, M. Liu, C. Martin, J. Elefteriades, and W. Sun. A machine learning approach to investigate the relationship between shape features and numerically predicted risk of ascending aortic aneurysm. Biomechanics and modeling in mechanobiology, 16, 04 2017. L. Liang, M. Liu, C. Martin, and W. Sun. A deep learning approach to estimate stress distribution: a fast and accurate surrogate of finite-element analysis. Journal of The Royal Society Interface, 15, 01 2018. R. Limet, N. Sakalihasan, and A. Albert. Determination of the expansion rate and incidence of rupture of abdominal aortic aneurysms. Journal of vascular surgery : official publication, the Society for Vascular Surgery and International Society for Cardiovascular Surgery, North American Chapter, 14:540–8, 11 1991. J. Long, E. Shelhamer, and T. Darrell. Fully convolutional networks for semantic segmenta- tion. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 3431–3440, 2015. Q. Long, X. Xu, M. Collins, T. Griffith, and M. Bourne. The combination of magnetic resonance angiography and computational fluid dynamics: A critical review. Critical reviews in biomedical engineering, 26:227–74, 02 1998. A. L. Marsden. Optimization in cardiovascular modeling. Annual Review of Fluid Mechanics, 46, 12 2013. A. L. Marsden, M. Wang, J. Dennis, and P. Moin. Suppression of vortex-shedding noise via derivative-free shape optimization. Physics of Fluids - PHYS FLUIDS, 16, 10 2004. A. L. Marsden, M. Wang, J. Dennis, and P. Moin. Trailing-edge noise reduction using derivative-free optimization and large-eddy simulation. Journal of Fluid Mechanics, 572: 13 – 36, 02 2007. A. L. Marsden, J. Feinstein, and C. Taylor. A computational framework for derivative-free optimization of cardiovascular geometries. Computer Methods in Applied Mechanics and Engineering, 197:1890–1905, 04 2008. G. Martufi, M. Auer, J. Roy, J. Swedenborg, N. Sakalihasan, G. Panuccio, and T. Gasser. Multidimensional growth measurements of abdominal aortic aneurysms. Journal of vas- cular surgery, 58, 04 2013. G. Matheron. Principles of geostatistics. Economic Geology, 58(8):1246–1266, 12 1963. Y. Matsuda, K. Takanashi, J. Takasu, N. Morooka, and Y. Inagaki. Expansion rate of thoracic aortic aneurysms and influencing factors. Chest, 102(2):461 – 466, 1992. R. Memisevic and G. E. Hinton. Unsupervised learning of image transformations. In 2007 IEEE Conference on Computer Vision and Pattern Recognition, pages 1–8. IEEE, 2007. 114 F. Meng, Z. Lu, M. Wang, H. Li, W. Jiang, and Q. Liu. Encoding source language with convolutional neural network for machine translation. In Proceedings of the 53rd An- nual Meeting of the Association for Computational Linguistics and the 7th International Joint Conference on Natural Language Processing (Volume 1: Long Papers), pages 20–30, Beijing, China, July 2015. Association for Computational Linguistics. S. J. Mousavi, S. Farzaneh, and S. Avril. Patient-specific predictions of aneurysm growth and remodeling in the ascending thoracic aorta using the homogenized constrained mixture model. Biomechanics and Modeling in Mechanobiology, 18:1895–1913, 2019. R. M. Neal. Priors for infinite networks. In Bayesian Learning for Neural Networks, Lecture Notes in Statistics, vol 118, Berlin, Heidelberg, 1996. Springer, New York, NY. S. A. Niederer, Y. Aboelkassem, C. D. Cantwell, C. Corrado, S. Coveney, E. M. Cherry, T. Delhaas, F. H. Fenton, A. V. Panfilov, P. Pathmanathan, G. Plank, M. Riabiz, C. H. Roney, R. W. dos Santos, and L. Wang. Creation and application of virtual patient cohorts of heart models. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 378(2173):20190558, 2020. R. A. Nishimura, C. M. Otto, R. O. Bonow, B. A. Carabello, J. P. Erwin, R. A. Guyton, P. T. O’Gara, C. E. Ruiz, N. J. Skubas, P. Sorajja, T. M. Sundt, and J. D. Thomas. 2014 aha/acc guideline for the management of patients with valvular heart disease. Journal of the American College of Cardiology, 63(22):e57–e185, 2014. D. Oliveira, J. Srinivasan, D. Espino, K. Buchan, D. Dawson, and D. Shepherd. Geometric description for the anatomy of the mitral valve: A review. Journal of Anatomy, 237(2): 209–224, 2020. C. M. Otto, R. A. Nishimura, R. O. Bonow, B. A. Carabello, J. P. Erwin, F. Gentile, H. Jneid, E. V. Krieger, M. Mack, C. McLeod, P. T. O’Gara, V. H. Rigolin, T. M. Sundt, A. Thompson, and C. Toly. 2020 acc/aha guideline for the management of patients with valvular heart disease: A report of the american college of cardiology/american heart association joint committee on clinical practice guidelines. Circulation, 143(5):e72–e227, 2021. K. I. Paraskevas, D. P. Mikhailidis, V. Andrikopoulos, N. Bessias, and S. P. R. F. Bell. Should the size threshold for elective abdominal aortic aneurysm repair be lowered in the endovascular era? yes. Angiology, 61:617–9, 2010. A. Parr, M. McCann, B. Bradshaw, A. Shahzad, P. Buttner, and J. Golledge. Thrombus volume is associated with cardiovascular events and aneurysm growth in patients who have abdominal aortic aneurysms. Journal of vascular surgery, 53:28–35, 10 2010. P. Perdikaris and G. Karniadakis. Model inversion via multi-fidelity bayesian optimization: A new paradigm for parameter estimation in haemodynamics, and beyond. Journal of The Royal Society Interface, 13:20151107, 05 2016. 115 P. Perdikaris, D. Venturi, J. Royset, and G. Karniadakis. Multi-fidelity modelling via recur- sive co-kriging and gaussian-markov random fields. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science, 471:20150018, 07 2015. T. Pham and W. Sun. Material properties of aged human mitral valve leaflets. Journal of biomedical materials research. Part A, 102, 08 2014. A. Quarteroni and G. Rozza. Optimal control and shape optimization of aorto-coronaric bypass anastomoses. Mathematical Models and Methods in Applied Sciences, 13, 11 2011. R. Rizzo, W. McCarthy, S. Dixit, M. Lilly, V. Shively, W. Flinn, and J. Yao. Collagen types and matrix protein content in human abdominal aortic aneurysms. Journal of vascular surgery : official publication, the Society for Vascular Surgery and International Society for Cardiovascular Surgery, North American Chapter, 10:365–73, 11 1989. J. Sacks, W. J. Welch, T. J. Mitchell, and H. P. Wynn. Design and analysis of computer experiments. Statistical Science, 4(4):409–423, 1989. S. Sankaran and A. L. Marsden. The impact of uncertainty on shape optimization of idealized bypass graft models in unsteady flow. Physics of Fluids, 22, 12 2010. S. Sankaran, J. D. Humphrey, and A. L. Marsden. An efficient framework for optimiza- tion and parameter sensitivity analysis in arterial growth and remodeling computations. Computer methods in applied mechanics and engineering, 256:200–210, 04 2013. S. Sarkar, S. Mondal, M. Joly, M. Lynch, S. Bopardikar, R. Acharya, and P. Perdikaris. Multifidelity and multiscale bayesian framework for high-dimensional engineering design and calibration. Journal of Mechanical Design, 141:1–24, 09 2019. C. Scotti, A. Shkolnik, S. Muluk, and E. Finol. Fluid-structure interaction in abdominal aortic aneurysms: Effects of asymmetry and wall thickness. Biomedical engineering online, 4:64, 02 2005. M. Sermesant, P. Moireau, O. Camara, J. Sainte-Marie, R. Andrian, R. Cimrman, D. Hill, D. Chapelle, and R. Razavi. Cardiac function estimation from mri using a heart model and data assimilation: Advances and difficulties. Medical image analysis, 10:642–56, 09 2006. S. Seyedsalehi, L. Zhang, J. Choi, and S. Baek. Prior distributions of material parameters for bayesian calibration of growth and remodeling computational model of abdominal aortic wall. Journal of biomechanical engineering, 137, 07 2015. A. Siika, M. L. Liljeqvist, R. Hultgren, T. Gasser, and J. Roy. Aaa rupture often occurs outside the maximal diameter region and is preceded by rapid local growth and an in- creased biomechanical rupture risk index. European Journal of Vascular and Endovascular Surgery, 50:390, 09 2015. R. Spilker and C. Taylor. Tuning multidomain hemodynamic simulations to match physio- logical measurements. Annals of biomedical engineering, 38:2635–48, 03 2010. 116 N. Srinivas, A. Krause, S. Kakade, and M. Seeger. Information-theoretic regret bounds for gaussian process optimization in the bandit setting. IEEE Transactions on Information Theory - TIT, 58:3250–3265, 05 2012. A. H. Stroud and D. Secrest. Gaussian quadrature formula. Mathematics of Computation, 21, 01 1966. M. J. Sweeting and S. G. Thompson. Making predictions from complex longitudinal data, with application to planning monitoring intervals in a national screening programme. Journal of the Royal Statistical Society: Series A (Statistics in Society), 175(2):569–586, 2012. J. M. Szafron, A. B. Ramachandra, C. K. Breuer, A. L. Marsden, and J. D. Humphrey. Optimization of tissue-engineered vascular graft design using computational modeling. Tissue Engineering Part C: Methods, 06 2019. C. Taylor, M. Draney, J. Ku, D. Parker, B. Steele, K. Wang, and C. Zarins. Predictive medicine: Computational techniques in therapeutic decision-making. Computer aided surgery : official journal of the International Society for Computer Aided Surgery, 4:231– 47, 01 1999. S. G. Thompson, H. A. Ashton, L. Gao, M. J. Buxton, and R. A. P. Scott. Final follow-up of the Multicentre Aneurysm Screening Study (MASS) randomized trial of abdominal aortic aneurysm screening. British Journal of Surgery, 99(12):1649–1656, 2012. K. A. Vardulaki, T. Prevost, N. Walker, N. Day, A. B. M. Wilmink, C. R. G. Quick, H. A. Ashton, and R. A. P. Scott. Growth rates and risk of rupture of abdominal aortic aneurysms. The British journal of surgery, 85:1674–80, 02 1999. M. Vega de Céniga, R. Gómez, L. Estallo, L. Rodrı́guez, M. Baquer, and A. Barba. Growth rate and associated factors in small abdominal aortic aneurysms. European Journal of Vascular and Endovascular Surgery, 31(3):231–236, 2006. F. Veronesi, C. Corsi, L. Sugeng, E. Caiani, L. Weinert, V. Mor-Avi, S. Cerutti, C. Lam- berti, and R. Lang. Quantification of mitral apparatus dynamics in functional and is- chemic mitral regurgitation using real-time 3-dimensional echocardiography. Journal of the American Society of Echocardiography : official publication of the American Society of Echocardiography, 21:347–54, 05 2008. J. Villadsen and M. L. Michelsen. Solution of differential equation models by polynomial approximation, volume 7. Prentice-Hall Englewood Cliffs, NJ, 1978. K. Y. Volokh and D. Vorp. A model of growth and rupture of abdominal aortic aneurysm. Journal of biomechanics, 41:1015–21, 02 2008. P. Watton and N. Hill. Evolving mechanical properties of a model of abdominal aortic aneurysm. Biomechanics and modeling in mechanobiology, 8:25–42, 01 2008. 117 C. K. I. Williams. Computing with infinite networks. In Proceedings of the 9th Interna- tional Conference on Neural Information Processing Systems, NIPS’96, pages 295–301, Cambridge, MA, USA, 1996. MIT Press. J. S. Wilson, S. Baek, and J. D. Humphrey. Importance of initial aortic properties on the evolving regional anisotropy, stiffness and wall thickness of human abdominal aortic aneurysms. Journal of the Royal Society, Interface / the Royal Society, 9:2047–58, 04 2012. J. S. Wilson, S. Baek, and J. D. Humphrey. Parametric study of effects of collagen turnover on the natural history of abdominal aortic aneurysms. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science, 469(2150), 2013. W. Yang, J. Feinstein, and A. L. Marsden. Constrained optimization of an idealized y-shaped baffle for the fontan surgery at rest and exercise. Computer Methods in Applied Mechanics and Engineering, 36):2135–2149, 07 2010. B. A. Zambrano, H. Gharahi, C. Lim, F. A. Jaberi, J. Choi, W. Lee, and S. Baek. Association of intraluminal thrombus, hemodynamic forces, and abdominal aortic aneurysm expansion using longitudinal CT images. Annals of Biomedical Engineering, 44:1502–14, 2015. B. A. Zambrano, N. McLean, X. Zhao, J. Tan, L. Zhong, C. Figueroa, L. C. Lee, and S. Baek. Image-based computational assessment of vascular wall mechanics and hemodynamics in pulmonary arterial hypertension patients. Journal of Biomechanics, 68, 12 2017. S. Zeinali-Davarani and S. Baek. Medical image-based simulation of abdominal aortic aneurysm growth. Mechanics Research Communications, 42:107–117, 06 2012. S. Zeinali-Davarani, A. Sheidaei, and S. Baek. A finite element model of stress-mediated vascular adaptation: Application to abdominal aortic aneurysms. Computer methods in biomechanics and biomedical engineering, 14:803–17, 04 2011. L. Zhang, Z. Jiang, J. Choi, C. Lim, T. Maiti, and S. Baek. Patient-specific prediction of abdominal aortic aneurysm expansion using bayesian calibration. IEEE Journal of Biomedical and Health Informatics, PP:1–1, 01 2019. Y. Zhou, Y. Wan, S. Roy, C. Taylor, C. R. Wanke, D. Ramamurthy, and J. Xie. Multivariate probabilistic collocation method for effective uncertainty evaluation with application to air traffic flow management. Systems, Man, and Cybernetics: Systems, IEEE Transactions on, 44(10):1347–1363, 2014. M. Zulliger. A strain energy function for arteries accounting for wall composition and struc- ture. Journal of Biomechanics, 37:989–1000, 07 2004. 118