Patient-Specific Modelling and Parameter Optimisation to Simulate Dilated Cardiomyopathy in Children

Purpose Lumped parameter modelling has been widely used to simulate cardiac function and physiological scenarios in cardiovascular research. Whereas several patient-specific lumped parameter models have been reported for adults, there is a limited number of studies aiming to simulate cardiac function in children. The aim of this study is to simulate patient-specific cardiovascular dynamics in children diagnosed with dilated cardiomyopathy, using a lumped parameter model. Methods Patient data including age, gender, heart rate, left and right ventricular end-systolic and end-diastolic volumes, cardiac output, systolic and diastolic aortic pressures were collected from 3 patients at Great Ormond Street Hospital for Children, London, UK. Ventricular geometrical data were additionally retrieved from cardiovascular magnetic resonance images. 23 parameters in the lumped parameter model were optimised to simulate systolic and diastolic pressures, end-systolic and end-diastolic volumes, cardiac output and left and right ventricular diameters in the patients using a direct search optimisation method. Results Difference between the haemodynamic parameters in the optimised cardiovascular system models and clinical data was less than 10%. Conclusion The simulation results show the potential of patient-specific lumped parameter modelling to simulate clinical cases. Modelling patient specific cardiac function and blood flow in the paediatric patients would allow us to evaluate a variety of physiological scenarios and treatment options.


INTRODUCTION
Lumped parameter modelling has been widely adopted to simulate cardiac function and physiological scenarios in cardiovascular research as they allow to study complex conditions in a simplified manner. For instance, heart failure has been simulated using this method by reducing systolic elastance, which describe the relationship between left ventricular pressure and volume. 2 Ventricular wall mechanics have been modelled using active and passive stresses on a fibre which in turn have been utilised to simulate and heart failure and mechanical circulatory support in adults. 5 Multiscale models which include cellular, protein and organ level dynamics have also been utilised to simulate heart failure in adults. 1 Atrial function in adults have been simulated using a multi-scale model. 25 Dynamics of heart valves have been modelled to simulate blood flow through the healthy and diseased heart valves in adults. 16 Lumped parameter models simulating heart failure or heart valve dynamics have been utilised to evaluate treatment techniques such heart pump support in adults. For instance, Gross et al. 12 evaluated the effect of left ventricular assist device speed increase during exercise in a model simulating heart failure. Liu et al. 20 evaluated varying speed modulation algorithms for continuous-flow left ventricular assist devices based on a cardiovascular coupling numerical model. Mod-elling of cardiac function and simulation of physiological cases in adults have been studied extensively and detailed information can be found in the literature. 18 Lumped parameter models simulating patientspecific cardiac dynamics in adults have already been developed. For instance, Ellwein et al. 7 developed a patient-specific model of cardiovascular and respiratory systems during hypercapnia. Pope et al. 26 showed how sensitivity analysis and subset selection can be employed in a cardiovascular model to estimate system parameters for healthy young and elderly subjects. Neal and Bassingthwaighte 22 developed a method to estimate subject-specific cardiac output and blood volume during haemorrhage. There is a limited number of studies aiming to simulate physiological cases in children where there is a need for such models. Schiavazzi et al. 29 estimated parameters of patient-specific in single-ventricle lumped circulation models under uncertainty. Goodwin et al. 10 developed a model simulating infant cardiovascular physiology for educational purposes. Sa´Couto et al. 28 developed a cardiovascular system model to simulate neonatal physiology. Giridharan et al. 9 developed a computer model of the paediatric circulatory system to test ventricular assist devices. Petukhov and Telyshev 24 developed a cardiovascular system model for paediatric patients with congenital heart diseases. Bozkurt 3 developed a cardiovascular system model capable of simulating cardiac function and heart chamber dimensions over a cardiac cycle in adults and children.
Mostly, average models describing generalised physiology in adults and children have been utilised to simulate physiological scenarios such as heart failure. However, causes of heart failure vary significantly from patient to patient and indeed, several different therapies, including digoxin, beta-blockers, diuretics or angiotensin-converting enzyme inhibitors, are used to improve heart failure. 14 Moreover, in children, clinical decisions are often dependant on results recorded in adult heart failure trials. 6 Therefore, treatment is negatively affected as there are substantial difficulties in adapting adult data to paediatric patients. Patientspecific numerical models simulating cardiac function may help to evaluate treatment options for children. Our previous model 3 presents a concept model in which the parameters were not optimized to simulate patient-specific clinical data. It simulates hemodynamic variables in adults and children within a physiological interval. The aim in this study is to simulate patient-specific cardiovascular dynamics in children diagnosed with dilated cardiomyopathy using a new optimisation strategy which can be applied without computational cost for parameter identification.

MATERIALS AND METHODS
The following demographics and hemodynamic parameters were collected from 3 patients diagnosed with dilated cardiomyopathy at Great Ormond Street Hospital for Children, London, UK, by searching the clinical database: age, gender, heart rate, cardiac output, left and right ventricular end-systolic and end-diastolic volumes, systolic and diastolic aortic pressures. The patients were selected to simulate mild and severe dilated cardiomyopathy in the left ventricle and biventricular dysfunction.
Cardiac magnetic resonance images (1.5-Tesla Magnetom Avanto scanner, Siemens Medical Solutions, Erlangen, Germany) were reviewed for these patients. The steady-state free precession sequences acquired at both end-diastole, and end-systole, were used to measure cardiac dimensions in Simpleware ScanIP 2018 (Synopsis, CA, USA). Four-chamber view was used to measure left and right ventricular long axis length and right ventricular basal diameter; the short-axis view was used to measure left ventricular lateral-septal diameter.

Patient Information
In children, dilated cardiomyopathy is diagnosed considering z-scores for ventricular ejection fraction, ventricular volume and ventricular end-diastolic and end-systolic diameters indexed on body surface area and measurements of the ventricular systolic function (z-score>2). 19,31 Patient 1 (male, age at scan = 8 years 7 months, body surface area = 0.94 m 2 ) was diagnosed with mild dilated cardiomyopathy. There were no signs of left ventricular hypertrophy, atrial dilation, atrio-ventricular valve regurgitation, and outflow tracts' obstruction. Both aortic and pulmonary valves were competent, with no stenosis. The right ventricle had normal global systolic function and volume. Cardiac output was 5.3 L/min, and systolic and diastolic pressures were 109 mmHg and 60 mmHg, respectively. End-diastolic volume was 106 mmHg (z-score=2.1). Patient 2 (female, age at scan = 10 years 8 months, body surface area = 0.94 m 2 ) was diagnosed with severe dilated cardiomyopathy. She presented with global thinning of the myocardium with no left ventricular hypertrophy. There was no left atrial dilatation and no mitral regurgitation on volume measurements. There was no left ventricular outflow tract obstruction and the aortic valve had normal function. The right atrium was not dilated. There was no tricuspid valve regurgitation. There was no right ventricular outflow tract obstruction. The pulmonary valve was competent with no stenosis. Cardiac output was 2.6 L/min, and systolic and diastolic pressures were 84 mmHg and 49 mmHg, respectively. Patient 3 (male, age at scan = 15 years 2 months, body surface area = 1.44 m 2 ) was diagnosed with severe cardiomyopathy with biventricular dysfunction and raised pulmonary pressure. The right atrium was dilated, whilst the left atrial chamber remained normal. There was 37% tricuspid valve regurgitation. The right ventricular myocardium was moderately hypertrophied, with globally poor contractility and evidence of severely impaired diastolic function. The right ventricle was dilated at both end-diastole and end-systole. There was no left ventricular hypertrophy. The left ventricular chamber was of low-normal body surface area indexed volume at end-diastole, and normal volume at end-systole, appearing under-filled, with a global impairment of systolic function. The pulmonary valve was unobstructed whilst there was 11% pulmonary valve regurgitation. The aortic valve was competent. Cardiac output was 2.6 L/min, and systolic and diastolic pressures were 118 mmHg and 64 mmHg, respectively.

Cardiovascular System Model
The cardiovascular system model used in this study simulates pressures, volumes and diameters in the heart chambers, flow rate through the heart valves, and pressures and flow rates in the systemic and pulmonary circulations. The left ventricular pressure (p lv ) was described using the active and passive contraction components (p lv,a , p lv,p ). The ventricular active pressure component was described using systolic ventricular elastance (E es,lv ), ventricular volume and zero-pressure volume (V lv , V lv,0 ), and the activation function (f act,lv ). The ventricular passive pressure component (p lv,p ) was modelled using an exponential relationship including volume (V lv ) and additional coefficients (A, B). Detailed information about the cardiovascular system model can be found in. 3 p lv ¼ p lv;a þ p lv;p ð1Þ p lv;a ¼ E es;lv ðV lv À V lv;0 Þf act;lv ðtÞ ð 2Þ Left ventricular volume (V lv ) was described using the left ventricular radius (r lv ), long axis length (l lv ) and an additional coefficient (K lv ), which includes effects of the contraction in the long axis and scales the proportion between the left ventricular radius and volume over a cardiac cycle. Change of the left ventricular radius (r lv ) were described utilising the flow rates through the aortic and mitral valves (Q av , Q mv ), the left ventricular volume (V lv ) and long axis length (l lv ), and the coefficient K lv .
The left atrial pressure (p la ) and volume (V la ) relationship was described using the elastance (E la ).
The left atrial volume (V la ) was described using the left atrial radius (r la ) and long axis length (l la ), and an additional coefficient (K la ). Changes of the left atrial radius (r la ) were described utilising the flow rates through the mitral valve and pulmonary vein (Q mv , Q vp ), the left atrial volume (V la ) and long axis length (l la ), and the coefficient K la .
The right atrium and ventricle were modelled in the same way as the left atrium and ventricle, using the different parameter values.
Heart valves were modelled using the pressure across the valve and the characteristic resistances for the forward and backward flows (R mv,f , R mv,b ) to simulate regurgitant valve flow in the patients. The flow rate expression through the mitral valve is provided below.
The other heart valves were modelled in a similar way using the appropriate parameter values.
The circulatory system included the aorta, the systemic arteries and veins, and the pulmonary arteries and veins. Blood flow in the circulatory system was described using a lumped parameter model, which included electrical analogues for resistance (R), compliance (C) and inertia (L). The aortic blood pressure (p ao ) and flow rate signals (Q ao ) are provided below.
here Q av and p as represent the aortic valve flow rate and systemic arteriolar pressure, respectively, and C ao , R ao and L ao representing compliance, resistance and inertance in the aorta, respectively. The electric-analogue of the cardiovascular system model is represented in Figure 1. Ejection fraction (EF) in the cardiovascular system models were calculated as given below.
Here, SV and EDV are the ventricular stroke and end-diastolic volumes.

Parameter Optimisation
Parameters were optimised for each patient to simulate systolic and diastolic pressures, end-systolic and end-diastolic volumes, cardiac output, and left and right ventricular diameters according to the retrieved clinical data, using a direct search method. Flowchart describing the optimisation algorithm is given in Figure 2.
The objective function (f) in the optimisation process aimed to minimise the difference in mean arterial pressure (MAP) and cardiac output (CO) between clinical data and the cardiovascular system model. The optimisation algorithm was designed to find a quick and effective solution for parameter estimation problem with multiple variables such as in the presented cardiovascular system model. The utilised optimisation method evaluates multiple variables in each simulation and effectively finds optimal solutions without computational cost.
The utilised electric analogue model of blood vessels simulates blood flow by solving ordinary differential equations for pressures and flow rates in each compartment. Heart chambers simulate pressure and volume which in turn affects the cardiac output the cardiovascular system. Error functions and objective function includes flow rate (CO) and pressure (MAP) terms to ensure that objective function is sensitive to the parameter set in the x vector. MAP in the cardiovascular system model is the time-averaged aortic pressure (p ao ) over a cardiac cycle in each patient.
here f MAP and f CO are the relative error functions for mean arterial pressure and cardiac output. The parameter set (x) included left and right ventricular end-systolic elastances (E es,lv , E es,rv ), left and right ventricular zero-pressure volumes (V 0,lv , V 0,rv ), the coefficients used in the passive pressure components of the ventricular pressures (A lv , A rv , B lv , B rv ), the compliances and resistances of the aorta, systemic arteries, and veins, pulmonary artery, arterioles and veins (R ao , C ao , R as , C as , R vs , C vs , R po , C po , R ap , C ap , R vp , C vp ) and the circulating blood volume (V blood ).
Upper and lower bounds of the optimised parameters (x ub , x lb ) were chosen around the values from the concept model published in. 3 Step-size (h) in the optimisation was defined using x ub , x lb and iteration number (n) in each simulation, and the parameters were updated using the initial lower bounds of the parameters and the step size, except for the coefficients A and B. These were updated starting with the upper bounds. Because arterial pressure and cardiac output decrease with the increasing A and B.
Step size (h) and parameter updates are given in the Equations 17-19.
The optimisation algorithm was set to complete 10 iterations in one simulation. The upper and lower bounds of the parameters were updated according to the values of the objective and error functions. The simulations were repeated to complete 10 iterations in each simulation with updated upper and lower bounds  The left and right ventricular diameters were estimated by optimising the coefficients in K lv and K rv .
The error function for the diastolic and systolic diameters of the left and right ventricles were defined as: Upper and lower bounds of the k vector (k ub , k lb ) were chosen around the values from the concept model published in. 3 Again, the step-size (h) in the optimisation was defined using k ub , k lb , and the iteration number (n) in each simulation and the parameters were updated using the initial lower bounds of the parameters and the step size. K lv and K rv do not affect any other variable but ventricular diameters, depending on the ventricular volumes. Therefore, first ventricular volumes were optimized and then K lv and K rv were tuned according to the ventricular volumes. Initial upper and lower bound of the parameters in the k vector are reported in Table 2.
Blood volume optimised in this study is the circulating blood volume. Total blood volume in children for the studied ages changes between 2000 mL and 4000 mL. 27 However, most of the blood volume consists of unstressed blood and constitutes a blood volume reserve. 11 Circulating blood volume is less than the total blood volume, therefore, optimised values for the circulating blood volume were less than the total blood volume in the cardiovascular system models.
Patient-specific heart rates and left and right ventricular long axis lengths (l v ) were used as input values TABLE 1. Initial upper and lower bounds of the parameters in the x and k vectors in the cardiovascular system models (E, V, A and B represent elastance, volume and the parameters used in the passive properties in the ventricle models, R, C and L represent resistance, compliance and inertance of the blood vessels, K is a coefficient in the ventricle models, subscripts es, lv, rv, 0 represent end-systole, left and right ventricles and initial value, ao, as, vs, po, ap, vp and blood represent aorta, systemic arteries and veins, pulmonary artery, arterioles and veins, and blood).  in the model. Ventricular long axis lengths were kept constant in the simulations and calculated using systolic and diastolic long axis lengths (l v,sys , l v,dias ) in each patient as: The parameters used in in the compartments describing atria, forward flow rate characteristics of the heart valves and inertance in the blood vessels were kept the same in all the models values because there were no data available to optimise them. Heart valve characteristics were assumed to be the same unless there were insufficiency. Moreover, inertial properties of blood had negligible effects on the CO and MAP. The parameter values that were not optimised for each specific case are reported in Table 2.
The optimisation and simulation processes were carried out in Matlab Simulink R2017a. The set of equations was solved using the ode15s solver. The maximum step size was 1e-3 s and the relative tolerance was set to 1e-3.

RESULTS
Haemodynamic variables collected from the clinical database and those from the numerical models simulating each patient specific cardiac function are provided in Table 3.
Haemodynamic values from the patients and numerical models were generally in good agreement. The largest differences were for blood pressures, in Patient 2 aortic pressure in diastole (6 mmHg); for ventricular volumes, in Patient 3 for the left ventricle in diastole (6 mL); for stroke volume, in Patient 1 for the right ventricle (7 mL). Ventricular diameters were also in good agreement, with the largest difference being 0.3 cm in diastolic right ventricular diameter in Patient 3. The upper and lower bounds of the parameters in the x and k vectors obtained from the initial optimisation and used in the second optimisation are given in Table 4.
The difference between the upper and lower bounds were narrowed down after the first optimisation with the initial upper and lower bounds of the parameters. The second optimisation allowed us to calculate optimal values satisfying conditions defined in the objective function and the error functions for the parameters in the x and k vectors, as given in the Ta-TABLE 3. Input variables (HR, l lv , l rv ) and the simulated haemodynamic variables in the patients and simulation results in the numerical models simulating patient specific cardiac function (HR, MAP and CO represent heart rate, mean arterial pressure and cardiac output, p, V, SV, EF, D and l represent pressure, volume, stroke volume, ejection fraction, diameter and length, subscripts ao, lv, rv, sys and dias represent aorta, left and right ventricles, systolic and diastolic phases respectively).   Table 6.
The objective function f was £0.050 for all models in the optimisation. Moreover, the error function remained between -0.1 and 0.1. The plots of left and right ventricular pressures, and aortic and pulmonary arterial pressures from the three models are shown in Figure 3.
Systolic and diastolic aortic pressures were around 106 mmHg and 57 mmHg correlating with the patient data in Patient 1. Systolic right ventricle and pulmonary arterial pressures were around 35 mmHg whereas diastolic pulmonary arterial pressure was around 20 mmHg. Reduced systolic left ventricular and aortic pressures Patient 2 were simulated as in the patient data. Right ventricular and pulmonary arterial pressures were above 30 mmHg at systole, whilst diastolic pulmonary arterial pressure was slightly above 10 mmHg. Systolic left ventricular and aortic pressures were around 118 mmHg in Patient 3. Both systolic and diastolic pressures were remarkably high in the right ventricle showing the effects of both systolic and diastolic dysfunctions in the right ventricle. Left and right ventricular volumes and basal diameters extracted from the models are plotted over 2 cardiac cycles in Figure 4.
The left ventricular volume was larger than the right ventricular volume during the entire cardiac cycle in both Patient 1 and 2. The right ventricle was remarkably dilated in Patient 3, whereas the left ventricular volume remained below 100 mL over the cardiac cycle. Moreover, the right ventricular stroke volume was larger than the left ventricular stoke volume in Patient 3 because of the reverse flow through the tricuspid and pulmonary valves. Left ventricular basal diameter was higher than the right ventricular basal diameter in Patient 1 and 2, whilst right ventricular basal diameter was relatively high in Patient 3.
The left and right ventricular pressure-volume loops from the models are given in Figure 5.
The left ventricular pressure-volume loop shifted to the right in Patient 1 and Patient 2 due to systolic dysfunction in the left ventricle. In Patient 3, the right ventricular pressure-volume loop shifted to the right, indicating impaired right ventricular function. TABLE 4. The upper and lower bounds of the parameters in the x and k vectors in the second optimisation (E, V, A and B represent elastance, volume and the parameters used in the passive properties in the ventricle models, R, C and L represent resistance, compliance and inertance of the blood vessels, K is a coefficient in the ventricle models, subscripts es, lv, rv, 0 represent end-systole, left and right ventricles and initial value, ao, as, vs, po, ap, vp and blood represent aorta, systemic arteries and veins, pulmonary artery, arterioles and veins, and blood).

DISCUSSION
In this study, the parameters of a cardiovascular system model which included the heart, and the systemic and pulmonary circulations, were optimised to simulate patient specific cardiac function in 3 paediatric patients diagnosed with dilated cardiomyopathy. In children, dilated cardiomyopathy is diagnosed considering ventricular end-diastolic and end-systolic diameters indexed on body surface area, and measurements of the ventricular systolic function. 19 Left ventricular systolic elastance (E es,lv ) in adults is around 2.5 mmHg/mL whereas it is reported to be above 3 mmHg/mL in children between 6 and 15 ages. 13,15 Relatively high systolic elastance values in Patient 1 and Patient 3 show the effect of a mildly affected left ventricular function; whilst for Patient 2, diagnosed with severe left ventricular cardiomyopathy, left ventricular systolic elastance (E es,lv ) was the lowest (Table 5). A typical value right ventricular systolic elastance (E es,rv ) is around 1 mmHg/mL in adults, 4 whilst 1.4 mmHg/mL has been used for children between 8-12 years age. 3 1.040 mmHg/mL right ventricular elastance (E es,rv ) shows impaired right ventricular systolic function in Patient 2 (Table 5). Moreover, 0.338 mmHg/mL right ventricular elastance (E es,rv ) in Patient 3 ( Table 5) shows severely impaired right ventricular systolic function, correlating with the clinical findings.
The parameters A and B were used to describe the passive ventricular pressure-volume relation. The same intervals resulting in the same optimal values were used for the left and right ventricles in the models simulating cardiac function in Patient 1 and 2, in these patients, only systolic dysfunction, which is related to active contraction behaviour in the ventricles, has been reported. Therefore, left and right ventricular passive contraction properties were assumed to be similar in these patients. Conversely, Patient 3 was diagnosed with diastolic ventricular dysfunction in both ventricles. Therefore, relatively large intervals for the coefficients A and B were selected during the optimisation. Relatively low A and B values in the right ventricle model of Patient 3 shows severely impaired diastolic function in Patient 3, along with the systolic dysfunction.
Systemic vascular resistance increases in heart failure as a response to neurohumoral pathways to maintain perfusion level. 17 Increased systemic arteriolar resistance (R as ) can be seen in Patient 2 and 3 as the cardiac output and mean arterial pressure in these patient were relatively low. Upper and lower bounds of the other parameters were selected considering the concept model and haemodynamic variables such as ventricular volumes and arterial pressures reported in the patient clinical data. For instance, the interval between the upper and lower bounds in the ventricular zero-pressure volumes (V 0 ) were selected higher for the relatively high ventricular volumes clinically reported. Moreover, the right ventricular zero-pressure volumes (V rv,0 ) were higher in the concept model reported in. 3 Therefore, ventricular zero-pressure volumes (V 0 ) were the same for left and right ventricles in the cardiovascular system model simulating Patient 1 where right ventricular zero-pressure volumes (V rv,0 ) were higher as TABLE 5. Optimal values of the parameters in the x and k vector (E, V, A and B represent elastance, volume and the parameters used in the passive properties in the ventricle models, R, C and L represent resistance, compliance and inertance of the blood vessels, K is a coefficient in the ventricle models, subscripts es, lv, rv, 0 represent endsystole, left and right ventricles and initial value, ao, as, vs, po, ap, vp and blood represent aorta, systemic arteries and veins, pulmonary artery, arterioles and veins, and blood). in the cardiovascular system model simulating Patients 2 and 3. The x vector contained 21 parameters and the k vector contained 2 parameters. The optimal solution depends on the selected upper and lower bounds (ub, lb), and on the step-size (h). It may be possible to obtain the global solution or local minima which provide better solutions by changing the initial upper and lower bounds, step size or interval for the error functions. However, it should be noted that the aim of the study was to simulate clinical conditions within the specified limits for the error functions (f MAP and f CO ) and show that it could be possible to simulate clinical values in the patients using lumped parameter models. Also, relatively small step-sizes may allow to find an optimal solution in one run, whilst larger step-sizes updates the upper and lower bounds of the parameters and require running the optimisation algorithm more than once. Therefore, an initial simulation would be required for relatively large upper and lower bound intervals to narrow down the intervals for lower and upper bounds of the parameters. A numerical model utilising relatively large intervals for upper and lower bounds would take more time to find an optimal solution. The used optimisation method allowed us to optimise a parameter set containing several parameters without additional computational cost, as it utilised the parameter values simultaneously to find the optimal solution. Therefore, it is also suitable for low configuration computers.
Ventricular passive pressure component uses left ventricular volume as its independent variable. Although the ventricular pressures may not be zero at the specified zero-pressure volumes (V 0 ) values, the utilised model simulate left ventricular pressures within the physiological limits. Utilising also V 0 in the ventricular passive pressure component will allow ventricular pressures to become zero at the specified V 0 values. 21 Patient-specific lumped parameter models simulating cardiac function in children can be used to evaluate several different scenarios. In this study, dilated cardiomyopathy was considered, but the framework presented could be easily adapted to simulate for example restrictive or hypertrophic cardiomyopathy. 32 Another application would be simulations of therapies, such as ventricular assist device support that are becoming a viable treatment option for end-stage heart failure paediatric patients, 30 albeit often causing complica- tions which may require additional intervention. 8 Computational techniques have been utilised widely to evaluate clinical scenarios in children. 23 In those cases, where flow and pressure data cannot be retrieved from conventional imaging modalities, patient specific lumped parameter models, as presented in this study, can be used to simulate patient specific boundary conditions for finite element models.