Generative adversarial networks for construction of virtual populations of mechanistic models: simulations to study Omecamtiv Mecarbil action

Biophysical models are increasingly used to gain mechanistic insights by fitting and reproducing experimental and clinical data. The inherent variability in the recorded datasets, however, presents a key challenge. In this study, we present a novel approach, which integrates mechanistic modeling and machine learning to analyze in vitro cardiac mechanics data and solve the inverse problem of model parameter inference. We designed a novel generative adversarial network (GAN) and employed it to construct virtual populations of cardiac ventricular myocyte models in order to study the action of Omecamtiv Mecarbil (OM), a positive cardiac inotrope. Populations of models were calibrated from mechanically unloaded myocyte shortening recordings obtained in experiments on rat myocytes in the presence and absence of OM. The GAN was able to infer model parameters while incorporating prior information about which model parameters OM targets. The generated populations of models reproduced variations in myocyte contraction recorded during in vitro experiments and provided improved understanding of OM’s mechanism of action. Inverse mapping of the experimental data using our approach suggests a novel action of OM, whereby it modifies interactions between myosin and tropomyosin proteins. To validate our approach, the inferred model parameters were used to replicate other in vitro experimental protocols, such as skinned preparations demonstrating an increase in calcium sensitivity and a decrease in the Hill coefficient of the force–calcium (F–Ca) curve under OM action. Our approach thereby facilitated the identification of the mechanistic underpinnings of experimental observations and the exploration of different hypotheses regarding variability in this complex biological system. Supplementary Information The online version of this articlecontains supplementary material available 10.1007/s10928-021-09787-4.

group. The equation for the fraction of groups G XB engaged in force generation is where f G and g G rates are both functions of the fractions of troponin complexes with bound calcium (A), 12 and g G is also a function of the mean XB strain (s). The fraction of bound XBs (XB G ) within a group is 13 described as where the rate f XB is a constant parameter, and g XB is a function of the XB strain. The total fraction 15 of XBs in the force-generating state is equal to G XB × XB G . The rates f G and g G were formulated as 16 functions of A, where 0 ≤ ζ ≤ 1 andf G are model parameters, g G is a function of XB strain, and A 50 and n A are parameters 18 that define cooperativity in the interaction between troponin, tropomyosin state, and crossbridge dynamics 19 in the XB-groups. The equation for the mean XB strain is, where SL is the sarcomere length. The rates g XB and g G in (2) and (3) are functions of s. 21 g XB (s) =ḡ XB × s mod (s) where x 0 is the distortion due to XB power-stroke, andḡ G ,ḡ XB , α, and η are parameters of the model. 22 The length dependence for the thick filament is described by a piecewise polynomial function where ∧ and ∨ are min and max binary operators, respectively, λ is the stretch ratio of a sarcomere (we 24 assume a sarcomere length of 1.9 µm at λ = 1), λ ms0 > λ ms1 are the slopes of the piecewise linear function, 25 and λ mn0 < λ mn1 are the nodes defining the function discontinuities. The active tension developed by the 26 muscle (T a ) is then calculated as where T a is the active tension; S a is the scaling factor for the tension; x 0 is the step size of the myosin power stroke. Cardiac troponin complexes modulate myofilament contraction in a Ca-dependent manner by regulating the availability of myosin binding sites on actin. This is captured as in [9], where A H and A L refer to the fraction of high and low calcium affinity troponin sites, respectively, that have 28 Ca bound to their regulatory binding sites; k onT is the rate constant for binding; [Ca] is the concentration of 29 Ca; k offHT and k offLT are the rate constants for unbinding from the high and low affinity sites, respectively.

30
Similar to the thick filament, a function for length dependence is defined for the thin filament as The equation for intracellular calcium transient that drives the contraction is The myofilament model parameters were optimized to reproduce different in vitro experimental tests, such 37 as the rate of force redevelopment (Ktr) test (Fig 1A), isosarcometric test (Fig 1B), unloaded shortening 38 test ( Fig 1C) and isometric test (Fig 1D), typically carried out to characterize myofilament properties as 39 shown in Fig 1. The default model parameters are listed in Table 1 (Fig 1A), as observed experimentally [2,13]. The model simulations also showed the typical cooperative 43 relationship between force and calcium concentrations (i.e., significant increase in maximum force for small 44 changes in calcium levels) (Fig 1B) upon stimulation at fixed sarcomere lengths and Ca concentration. The 45 simulations also captured the experimentally observed [2] increases in maximal steady force and calcium 46 sensitivity at higher sarcomere lengths. Finally, the model simulation produced increases in isometric force 47 at larger sarcomere lengths ( Fig 1D).   The GAN architecture shown in Fig. 3 was used to train 3 generators to sample model parameters coherent 64 with the experimental data in both control and OM conditions. To stabilize training and help prevent 65 mode collapse, we also incorporated a reconstruction network that takes generated parameter samples 66 concatenated from all 3 generators as input, and aims to reproduce the 3 sets of base variables z. This 67 approach is similar to that used in VEEGAN [12], but simplified, as it includes only the 2 component of 68 the reconstructor loss, excludes the cross-entropy term and removes the dependence of the discriminator 69 on z.
Rec. loss Rec. net Fig. 3 Generative network for model parameter inference. The generator network is trained to transform random variables Z 1 , Z 2 , and Z 3 with base Gaussian distributions to random variables with densities q Xc,g (xc,g) and q X d ,g (x d,g ) as approximations of q Xc (xc) and q X d (x d ). The generator factorizes density by using 3 networks G 1 , G 2 , and G 3 . The network G 1 is responsible for parameters x 1 that do not change under the drug action. G 2 is responsible for parameters that are affected by the drug x 2,c and generates their values for the control group. G 3 is the same as G 2 , but for the group under action of the drug. The conditional dependence of x 2,d and x 2,c on x 1 is implemented by the input of samples from the base distribution Z 1 for G 1 to both G 2 and G 3 . Parameters are pushed through the model y = M (x) to obtain q Yc,g (yc,g) and q Y d ,g (y d,g ) as approximations of q Yc (yc) and q Y d (y d ). Discriminators D 1 , D 2 separates samples xc,g and x d,g from samples of the prior distribution of the parameters (uniform in our case). D 3 and D 4 are discriminators for model outputs. The additional structure of reconstruction network is added for GAN stabilization.

GAN loss functions
is maximized. Discriminators D1 and D2 (termed D X ) distinguish between samples from the prior over 74 mechanistic parameters X c,prior , X d,prior and samples generated by G1, G2 and G3, for which the standard is maximized. The reconstruction network R aims to reproduce the original base distribution Z from 77 samples generated by G, for which the squared loss is calculated The generator network G generates mechanistic parameter sets from the base variable Z, for which 79 losses are calculated from all D Y and D X according to The total loss for G is then the weighted sum loss which is minimized, where w Y = 1.0, w X = 0.1, and w R = 1.0 are used as default weights. 82 We used the Adam optimizer with step size of 0.0001 for G and R, and 0.00002 for D X , and D Y .

83
The β 1 and β 2 parameters of the Adam optimizer were set to default values of 0.9 and 0.999, respectively, 84 as suggested in [5]. Mini-batch size was 100. Training was performed in two stages. First, G, R and D X  Table 3. out 20% of the dataset as test data, and trained the neural network for 100 epochs. Fig. 4 shows that the 100 prediction performance of a trained surrogate on the test data was good, with the surrogate accurately 101 approximating the mechanistic model outputs. To test whether the results produced by the trained GAN architecture were robust given different trained 104 surrogates or initial conditions, we performed 5 trials of GAN training, once using a different surrogate model 105 across each trial, and once with a different random seed across each trial. The surrogates were trained using 106 the same network structure, training data, and method, but will exhibit differences in network weights 107 after training, and therefore some differences in predicted feature values, which will influence the GAN 108 training. Fig. 5 demonstrates that the results we obtained were robust to deviations caused by the use 109 of different surrogate models. Similarly, different initial conditions will influence the final trained weights 110 of the GAN networks, even with the same surrogate model, so we tested different initial conditions using 111 different random seeds. Fig. 6 shows that the parameters sampled by the trained GAN were also robust to 112 differences in initial conditions. 113 Fig. 5 Effect of multiple surrogates on parameters sampled by the GAN. Each 'trial' used a different surrogate model, trained using the same method and training data. Similar marginal distributions of parameters were sampled across trials for both control and OM conditions, indicating that the GAN training was robust to changes in the surrogate model. Fig. 6 Effect of multiple different random seeds on parameters sampled by the GAN. Each 'trial' used an identical surrogate model, but a different random seed for initialization of weights in the GAN. Similar marginal distributions of parameters were sampled across trials for both control and OM conditions, indicating that the GAN training was robust across trials with different initial conditions.

91
In previous work, we demonstrated an efficient method of training a low order model of ventricular mechanics 115 via congruency training to reproduce the global behavior of a 3-D finite element model by scaling the 116 outputs of just a single modeled cell [1]. Briefly, the low order models comprise the myofilament model 117 coupled with a transformation module whose parameters are sized by machine learning to match the low 118 order model outputs to features of the detailed finite element models. We demonstrated that simple linear 119 transformations between sarcomere strain (tension) and ventricular volume (pressure) were sufficient to 120 reproduce global pressure-volume outputs of 3-D finite element models. Here, we demonstrate our first 121 results from simulations of OM with the low order model of ventricular mechanics. Although the results 122 are preliminary, we include them here to outline a pipeline for translation of in vitro experiments into 123 predictions of whole organ behavior.

Models of left ventricle 125
In our pipeline, we translated the in vitro simulation results described in this paper into the hemodynamics 126 of the left ventricle (LV), simulating the effect of OM on the population of healthy (N) and heart failure 127 (HF) hearts. For these purposes, geometric information based on echocardiographic data was extracted from 128 retrospective medical records from healthy subjects (n=97) and patients affected by heart failure (n=261). 129 As in [3], left ventricular geometry was parameterized with 6 parameters: the outer radius at base 130 Rb, the length of the longitudinal semi-axis of the outer spheroid Z, the ventricular wall thicknesses at 131 base L and apex H, the sphericity/conicity of the spheroid e ∈ [0, 1], and the truncation angle Ψ 0 . In the 132 first step, we used a regression model to obtain 6 parameters that describe LV geometry based on several 133 echocardiographic features: end-diastolic volume (EDV), end-diastolic diameter, septal wall thickness, and 134 lateral wall thickness. The publicly available SunnyBrook CMR dataset [7] was employed to train the linear 135 regression: 136 1. The parameters Rb and Z were reformulated by subtracting L and H, respectively, in order to decouple 137 these variables and minimize their dependency, 138 2. All geometry was standardized in volume (100mL) via similarity transformations to remove ventricle 139 size variability. This standard volume was obtained by multiplying the parameters L, Rb, Z and H by 140 the cubic root of the ratio between the target and the actual volume. The dimensionless parameters Ψ 0 141 and e were preserved.

147
The LV model was coupled with a Windkessel model of blood circulation, as in [1]. We fit the LV model 148 parameters S a and P at to ensure an optimal match with the EDV and end-systolic volume (ESV). S a 149 and P at represent the overall efficiency of contraction and a simplified representation of atrial pressure, 150 respectively. To adapt parameters of the Windkessel model and achieve a central aortic blood pressure 151 within the normal range of 100-120 mmHg, the hemodynamic parameters of the 3-element Windkessel 152 model, i.e., aortic valve resistance R 1 , aortic resistance R 2 , distal resistance R 3 and capacitance C, were 153 scaled proportionally to the stroke volume of each patient. The optimization problem was solved with the 154 L-BFGS-B algorithm. All other parameters were taken from [1]. Average mismatches to the target variables 155 were smaller than 2%.

156
To reproduce the OM effect, we changed the parameters n A and A 50 as shown in Table 4. Comparison 157 analyses were carried out with respect to the pressure and volume traces, computed for the optimal set of 158 parameters of each patient.  Fig. 7 shows simulated global outputs after the model parameters were optimized to match ESV and EDV 161 for the N (blue color) and HF (orange color) cases. As expected, optimized values of the S a parameters 162 were significantly lower (mean -46%) for the HF compared to N. Optimized values of P at parameters did 163 not show differences (mean -5%) for HF compared to N. As shown in Fig. 7, the groups are divided by 164 indices and correspond to the observed features in HF and N cases.

159
165 Fig. 8 shows traces demonstrating the contractile effect of simulated OM (dashed lines) for a low (left) 166 and high (right) concentration of OM, compared to representative simulated HF traces (solid lines). As 167 observed, OM decreased ESV and EDV but did not affect dP/dt. This results shows that OM's inotropic 168 effect is limited by the impairment in diastolic function. 169 Fig. 9 shows changes from baseline to low (blue) and high (orange) concentrations of OM. The top 170 row shows that the main effect of OM on left ventricular function is associated with a decrease in EDV 171 and ESV. The difference between the HF and N groups is related to the size of the ventricles. Changes in 172 ejection fraction between groups are the same, and changes in stroke volume are greater for HF compared 173 to N (middle row). The bottom row shows that OM increases the ejection time for HF more than for N.