Neural network emulation of the human ventricular cardiomyocyte action potential for more efficient computations in pharmacological studies
Abstract
Computer models of the human ventricular cardiomyocyte action potential (AP) have reached a level of detail and maturity that has led to an increasing number of applications in the pharmaceutical sector. However, interfacing the models with experimental data can become a significant computational burden. To mitigate the computational burden, the present study introduces a neural network (NN) that emulates the AP for given maximum conductances of selected ion channels, pumps, and exchangers. Its applicability in pharmacological studies was tested on synthetic and experimental data. The NN emulator potentially enables massive speedups compared to regular simulations and the forward problem (find drugged AP for pharmacological parameters defined as scaling factors of control maximum conductances) on synthetic data could be solved with average rootmeansquare errors (RMSE) of 0.47 mV in normal APs and of 14.5 mV in abnormal APs exhibiting early afterdepolarizations (72.5% of the emulated APs were alining with the abnormality, and the substantial majority of the remaining APs demonstrated pronounced proximity). This demonstrates not only very fast and mostly very accurate AP emulations but also the capability of accounting for discontinuities, a major advantage over existing emulation strategies. Furthermore, the inverse problem (find pharmacological parameters for control and drugged APs through optimization) on synthetic data could be solved with high accuracy shown by a maximum RMSE of 0.22 in the estimated pharmacological parameters. However, notable mismatches were observed between pharmacological parameters estimated from experimental data and distributions obtained from the Comprehensive in vitro Proarrhythmia Assay initiative. This reveals larger inaccuracies which can be attributed particularly to the fact that small tissue preparations were studied while the emulator was trained on single cardiomyocyte data. Overall, our study highlights the potential of NN emulators as powerful tool for an increased efficiency in future quantitative systems pharmacology studies.
eLife assessment
This valuable prospective study develops a new tool to accelerate pharmacological studies by using neural networks to emulate the human ventricular cardiomyocyte action potential. The evidence supporting the conclusions is convincing, based on using a large and highquality dataset to train the neural network emulator. There are nevertheless a few areas in which the article may be improved through validating the neural network emulators against extensive experimental data. In addition, the article may be improved through delineating the exact speedup achieved and the scope for acceleration. The work will be of broad interest to scientists working in cardiac simulation and quantitative system pharmacology.
https://doi.org/10.7554/eLife.91911.3.sa0Introduction
Computer models of human physiology are becoming increasingly detailed and mature and the area of ventricular cardiomyocyte electrophysiology (EP) is one of the most advanced. The most updated models include fine representations of ion movements through various important channels, pumps, and exchangers, and take the complex handling of intracellular calcium accurately into account (Grandi et al., 2010; O’Hara et al., 2011; Tomek et al., 2019; Bartolucci et al., 2020). While these models have individual strengths and limitations in replicating different aspects of physiology, pathology, and pharmacology (Corrado et al., 2021; Amuzescu et al., 2021), their degree of credibility has reached a level that has led to an increasing number of applications in academia and beyond. This holds in particular for the pharmaceutical sector, where much effort is spent on using computer modeling to reduce traditional preclinical and clinical methodologies for assessing the efficacy and safety of novel drug candidates (Mirams et al., 2011; Passini et al., 2017; Li et al., 2019; Passini et al., 2021). To improve the regulatory assessment of a drug’s proarrhythmic potential, the Comprehensive in Vitro Proarrhythmia Assay (CiPA) was proposed in 2013 following a workshop at the US Food and Drug Administration (Sager et al., 2014; Colatsky et al., 2016; Strauss et al., 2021). A central component is a computer model of human ventricular cardiomyocyte EP that is coupled to a pharmacological model describing the interaction between a given drug and multiple arrhythmiarelevant channels (Dutta et al., 2017; Li et al., 2017; Li et al., 2019). For a given drug, experimental channel block data are collected to inform the pharmacological model and corresponding simulations of the action potentials (AP; time course of the transmembrane potential) are performed to predict the proarrhythmic risk based on a mechanistically motivated biomarker (Chang et al., 2017; Li et al., 2017; Li et al., 2019). The prediction is then compared with findings in experimental (Blinova et al., 2018) and clinical (Vicente et al., 2018) studies. To compute the drugged AP for given pharmacological parameters is a forward problem, while the corresponding inverse problem is to find pharmacological parameters for given control (before drug administration) and drugged AP. Some relevant examples for the latter have been presented by Bottino et al., 2006 who estimated pharmacological parameters from APs of canine Purkinje fibers and by Tveito et al. who estimated pharmacological parameters from AP biomarkers measured in human induced pluripotent stem cellderived cardiomyocytes (Tveito et al., 2018) and several animal ventricular cardiomyocytes (Tveito et al., 2020). Furthermore, Jaeger et al., 2021 identified the optimal polypharmacological treatment for recovering APs of mutant ventricular cardiomyocytes based on biomarkers of simulated wild type and mutant APs.
When the models are interfaced with experimental data, attention should be paid to the inherent uncertainty in the data that results from beattobeat variability (intrinsic variability), celltocell variability (extrinsic variability), and measurement errors (observational uncertainty) (Mirams et al., 2016). Uncertainty propagates through the given problem from APs to estimated parameters or from parameters to predicted APs and must be properly quantified to draw reliable conclusions from the results. Multiple methodologies exist for this purpose (Oakley and O’Hagan, 2004; Mirams et al., 2016; Sher et al., 2022) but usually require many simulations, which even for ordinary differential equation (ODE)based models of cardiomyocyte EP can become a significant computational burden when considering that each simulation includes a substantial number of beats to reach the model’s limit cycle, (also often called steady state). To overcome this problem, surrogate models have emerged which approximate (emulate) chosen outputs for given inputs multiple orders of magnitude faster. In line with uncertainty quantification literature, the cardiomyocyte EP model is from now on termed the ‘simulator’, whereas the surrogate model is termed the ‘emulator’. Earlier work has reported on an emulator based on linear interpolation of a multidimensional lookup table Mirams et al., 2014 and more recently, Gaussian process (GP) emulators have become popular. Their key advantage is that in and outputs are modeled as random distributions which allows for rapid sampling of the posterior distributions (Chang et al., 2015; Johnstone et al., 2016; Coveney and Clayton, 2018; Ghosh, 2018; Rasmussen, 2019; Coveney et al., 2021) and while outputs of recently published GP emulators were relevant biomarkers of the AP (Chang et al., 2015; Johnstone et al., 2016; Coveney and Clayton, 2018; Ghosh, 2018; Coveney et al., 2021), the emulation of the entire AP can also be realized, for example through dimensionality reduction techniques such as the principal component analysis or regressing statetransition models (Mohammadi et al., 2019). However, GP emulators are not well suited to capture discontinuities of the response surface as standard GP emulators assume a smooth and continuous response to changes in parameter values. Applying GPs for modeling discontinuous functions therefore remains a largely open problem. Thus, they may fail to capture AP abnormalities, which is a particular drawback for pharmacological studies where bifurcations in behavior such as early afterdepolarizations (EAD) can occur (Ghosh, 2018). To address this, Ghosh, 2018 presented a twostep approach for the emulation of the AP duration at 90% repolarization that first sets up a GP for the location of discontinuities and then fits separate GP emulators for the output of interest either side of these boundaries. In contrast, it has been proven that neural networks (NN) can approximate even discontinuous functions with arbitrary precision in theory (Hornik et al., 1989), while recent works using NNs show empirically promising results for modeling partial differential equations containing discontinuities (Jagtap et al., 2020). These features render NN emulators suitable emulation candidates and while Lei and Mirams, 2021 have recently investigated NN emulation of hERG channel kinetics, Jeong et al., 2023 proposed a neural network using AP shapes as input for the prediction of a drug’s proarrhythmic risk. However, to the best of our knowledge, NN emulators have not yet been used as surrogate for cardiomyocyte EP models.
The present study introduces NN emulation of the human ventricular cardiomyocyte AP and investigates the applicability in pharmacological studies. To this end, a NN emulator was developed based on data generated using a stateoftheart simulator (Tomek et al., 2019; Tomek et al., 2020) and the evaluation was done for forward and inverse problems on synthetic and experimental data.
Materials and methods
The methodology of this study including the development of the emulator and the evaluation is outlined in Figure 1.
Simulator
Request a detailed protocolThe simulator of Tomek et al., 2020 (ToRORddynCl simulator) was used. This is available in CARPentry (Vigmond et al., 2008) and was implemented based on the published CellML file for the endocardial subtype (ToRORd_dynCl_endo.cellml, https://github.com/jtmff/torord/tree/master/cellml; jtmff, 2020). Simulations were performed in CARPentry with the single cell tool bench. To compute the gating variables, the RushLarsen Method (Rush and Larsen, 1978) was employed, which uses an analytical solution assuming fixed voltage over a small timestep, and the remaining variables were computed by the Forward Euler method. To ensure low computational cost, we found the maximum solver and sampling time steps that still produce accurate results as follows. Various solver and sampling time steps were applied to generate APs and the biomarkers (AP biomarkers and abnormalities) used in this study were computed and compared with those that correspond to the minimum time steps (solver: $0.005\text{}\mathrm{m}\mathrm{s}$; sampling: $0.01\text{}\mathrm{m}\mathrm{s}$). We considered the time steps with only 2% relative difference for all AP biomarkers (solver: $0.01\text{}\mathrm{m}\mathrm{s}$; sampling: $0.05\text{}\mathrm{m}\mathrm{s}$) to offer a sufficiently good approximation. APs were stimulated at a pacing cycle length of $1000\text{}\mathrm{m}\mathrm{s}$ using $1\text{}\mathrm{m}\mathrm{s}$ long rectangular transmembrane current density pulses of $53\text{}\frac{\mu \mathrm{A}}{{\mathrm{c}\mathrm{m}}^{2}}$ (Tomek et al., 2019). To approach the simulator’s steady state, a series of 1000 stimuli were applied for each new parameter set starting from the initial states specified in the CellML file (when 1000 additional stimuli were applied, the maximum intracellular [Ca^{2+}], [Cl^{}], and [Na^{+}] changed by 1.5%, 0.15%, and 1.7%, respectively). Note, that the simulations can also be performed using opensource software such as (Clerx et al., 2016), OpenCOR (Garny and Hunter, 2015) and openCARP (Plank et al., 2021).
Emulator
Here, we present the architecture and capabilities of the implemented and trained emulator. Details on the practical implementation and a link to the public code can be found in Code & data availability.
Data
We generated three data sets in the study to train, validate and test the emulator performance.
Training/validation data (#1)
Request a detailed protocolThe first data set is the the supervised training data set, containing pairs of maximum conductance samples $\mathbf{x}$ and the corresponding AP ${\hat{V}}_{\mathrm{m}}(t)$ that was obtained from the simulator. Sobol’ sequences were used to generate 40,000 maximum conductance samples, containing 20,000 maximum conductance samples between 0% and 200% of the original values and 20,000 maximum conductance samples between 50% and 150% of the original values. The first covers a range that was considered plausible in terms of physiology and pathology (Britton et al., 2017; Tomek et al., 2019), and in terms of pharmacology (where full block is plausible). The latter covers a range that was considered particularly relevant in line with experimental calibration results presented in Tomek et al., 2019 and patch clamp measurements of channels that were exposed to 30 clinical drugs blocks in up to the fourfold of the maximum free therapeutic concentration were analyzed in agreement with the CiPA paradigm (Li et al., 2019; Crumb et al., 2016). The SALibSensitivity Analysis Library (Herman and Usher, 2017) was used in the entire study for the generation of samples based on Sobol’ sequences. For each maximum conductance sample, simulations were performed to obtain the corresponding 40,000 APs. APs with a transmembrane potential difference of more than 10% of the amplitude between $t=0$ and $1000\text{}\mathrm{m}\mathrm{s}$ (indicative of an AP that is far away from full repolarization) were excluded. This, however, applied to only 116 APs.
Starting from the original APs in data set #1, the data were first extended by $10\text{}\mathrm{m}\mathrm{s}$ from $t\in [0,\text{}1000]\text{}\mathrm{m}\mathrm{s}$ to $t\in [10,\text{}1000]\text{}\mathrm{m}\mathrm{s}$ to enable some extrapolation of ${V}_{\mathrm{m}}$ and hence a better alignment of the depolarization; for $t\in [10,\text{}0]\text{}\mathrm{m}\mathrm{s}$ the initial resting membrane potential ${V}_{\mathrm{m}}(0)$ was held constant. Then, the data were nonuniformly resampled from the original uniformly simulated APs, to emphasize the depolarization slope with a high accuracy while lowering the number of repolarization samples. For this purpose, we resamled the APs to $4\text{}\mathrm{k}\mathrm{H}\mathrm{z}$ for $t\in [20,\text{}5)\text{}\mathrm{m}\mathrm{s}$ (resting phase) and $10\text{}\mathrm{k}\mathrm{H}\mathrm{z}$ for $t\in [5,\text{}20)\text{}\mathrm{m}\mathrm{s}$ (depolarization phase) to $4\text{}\mathrm{k}\mathrm{H}\mathrm{z}$. The repolarization phase ($t\in [20,\text{}1000]\text{}\mathrm{m}\mathrm{s}$) was also resampled to $4\text{}\mathrm{k}\mathrm{H}\mathrm{z}$.
From the initial training data set, 20% were randomly chosen to be used for validation leaving 31908 pairs of maximum conductances and corresponding APs for training. Figure 2 shows processed APs that were used for training, validation and the APs excluded due to missing full repolarization as described above.
Synthetic test data (#2/#3)
Request a detailed protocolTwo sets of synthetic data were created using the simulator and each of the sets consisted of control and drug data with pairs of maximum conductances and corresponding APs.
The control data were the same in both sets. They were created using an experimentally calibrated population of 100 synthetic cardiomyocytes (Britton et al., 2013; Muszkiewicz et al., 2016; Gemmell et al., 2016). and to this end, Sobol’ sequences were used to generate samples of maximum conductances with values between 50% and 150% of the original values. Maximum conductance samples that produced APs without abnormalities (checked for the last two consecutive APs; see AP biomarkers and abnormalities), and with seven biomarker values (derived from the last AP; see AP biomarkers and abnormalities) in agreement with experimental ranges (Table 1) were included in the population. Please note that the experimental ranges were not derived from the data set described in Experimental data (#4).
Data set #2: The motivation for creating data set #2 was to evaluate the emulator on data of normal APs. Drug data were created using 100 synthetic drugs represented by a set of pharmacological parameters. Each synthetic drug was built to have four different targets, with all channels, pumps, and exchangers related to the emulator inputs considered as potential targets. To this end, 100 samples of four pharmacological parameters, each with values between 0.5 (50% block) and 1.5 (50% enhancement) were randomly generated. The synthetic drugs were applied to the entire synthetic cardiomyocyte population by scaling each of the relevant control maximum conductances with the corresponding pharmacological parameter. The samples that produced APs without abnormalities (checked for the last two consecutive APs; see AP biomarkers and abnormalities) were included in the data set. No sample was excluded and thus, the data set consisted of 100 control data pairs and 10,000 drug data pairs.
Data set #3: The motivation for creating data set #3 was to test the emulator on data of abnormal APs showing the repolarization abnormality EAD. This is considered a particularly relevant AP abnormality in pharmacological studies because of their role in the genesis of druginduced ventricular arrhythmia’s (Weiss et al., 2010). Drug data were created using 10 synthetic drugs with the hERG channel and the Cav1.2 channel as targets. To this end, 10 samples with pharmacological parameters for $\mathrm{G}}_{\mathrm{K}\mathrm{r}$ and $\mathrm{P}}_{\mathrm{C}\mathrm{a}$ (Table 2) were generated and the synthetic drugs were applied to the entire synthetic cardiomyocyte population by scaling $\mathrm{G}}_{\mathrm{K}\mathrm{r}$ and $\mathrm{P}}_{\mathrm{C}\mathrm{a}$ with the corresponding pharmacological parameter. Of the 1000 APs simulated, we discarded APs with a transmembrane potential difference of more than 10% of the amplitude between $t=0$ and $1000\text{}\mathrm{m}\mathrm{s}$ (checked for the last AP), indicative of an AP that is far away from fully repolarizing within $1000\text{}\mathrm{m}\mathrm{s}$. This left us with 950 APs, 171 of which exhibit EAD (see EAD classification).
Experimental data (#4)
Request a detailed protocolIn the experimental data set, APs were recorded in small right ventricular trabeculae and papillary tissue preparations that were isolated from healthy human hearts (Orvos, 2019). The hearts were obtained from organ donors whose hearts were explanted to obtain pulmonary and aortic valves for transplant surgery. Before cardiac explantation, organ donors did not receive medication apart from dobutamine, furosemide, and plasma expanders. Proper consent was obtained for use of each individual’s tissue for experimentation. The conventional microelectrode technique was used for AP recordings and all measurements were carried out at 37°C. Stimulation of APs was done at a pacing cycle length of $1000\text{}\mathrm{m}\mathrm{s}$ using a pair of platinum electrodes that provided rectangular current pulses of $2\text{}\mathrm{m}\mathrm{s}$ duration. To allow the preparations to equilibrate, stimuli were delivered for at least $60\text{}\mathrm{m}\mathrm{i}\mathrm{n}$ before the measurements started. Measurements were performed under control conditions and after administration of five channelblocking drugs at one concentration in multiple preparations. Drugs were cisapride ($30\text{}\mathrm{n}\mathrm{M}$), dofetilide ($10\text{}\mathrm{n}\mathrm{M}$), sotalol ($30\text{}\mu \mathrm{m}$), terfenadine ($1\text{}\mu \mathrm{m}$), and verapamil ($300\text{}\mathrm{n}\mathrm{M}$). The last 10 consecutive APs of each measurement were analyzed to quantify the beattobeat variability. Overall, the beattobeat variability was found to be small (standard deviation in all APs below $7\text{}\mathrm{m}\mathrm{V}$ before the peak due to time alignment mismatch and below $2\text{}\mathrm{m}\mathrm{V}$ after the peak) and thus, the last 10 consecutive APs of each measurement were averaged. In most of the preparations, the standard deviation between beats did not vary over time and thus, no temporal correlation of noise was assumed. Averaging also reduced the noise level. The experimental data set contained one pair of averaged control and drugged AP per preparation per drug. Pairs were excluded if the biomarker values (see AP biomarkers and abnormalities) of the control or the drugged AP were not in the range found in the training data (see Training). This applied to seven pairs and the final data set contained three pairs for cisapride, dofetilide, sotalol, and terfenadine, and one pair for verapamil. All measurements were performed at the University of Szeged, Hungary, and conformed to the principles of the Declaration of Helsinki. Experimental protocols were approved by the University of Szeged and by the National Scientific and Research Ethical Review Boards (No. 5157/1997 OEj and 49910/20101018EKU [339/PI/010]).
An overview of all utilized data is given in Table 3.
Architecture
Request a detailed protocolThe emulator takes maximum conductances of channels, pumps, and exchangers as inputs and computes the corresponding AP (${V}_{\mathrm{m}}(t)$) after the last stimulus as output. It was trained to represent human ventricular cardiomyocytes under control and drugged conditions and the inputs were selected based on two assumptions: (1) The kinetics of channels are preserved, while the number of channels, pumps, and exchangers vary due to different expression levels (Syed et al., 2005; Groenendaal et al., 2015; KroghMadsen et al., 2016). These numbers are captured in the simulator by the maximum conductance parameters (or permeability parameters but maximum conductance is used as general term here for the sake of simplicity) that determine the respective current densities; (2) Channels, pumps, and exchangers are potential (intended and unintended) drug targets and the interaction between drugs and their targets can be described by a scaling of the related maximum conductances (Brennan et al., 2009). The corresponding scaling factors are pharmacological parameters. These assumptions allowed us to focus on maximum conductances and we considered those as inputs which either the AP is sensitive to (G_{Na}, G_{NaL}, P_{Ca}, G_{to}, G_{K}r, G_{K1}, G_{NCX}, P_{NaK}) or which are related to common drug targets (G_{Na}, G_{NaL}, P_{Ca}, G_{to}, G_{Kr}, G_{Ks}, G_{K1}; Crumb et al., 2016) leading to the following selection: G_{Na}, G_{NaL}, P_{Ca}, G_{to}, G_{Kr}, G_{Ks}, G_{K1}, G_{NCX}, and P_{NaK}. AP sensitivity was quantified using a global sensitivity analysis (GSA; see.Appendix 2) and the inclusion threshold was a totaleffect Sobol’ sensitivity index (ST) above 0.1 with respect to any of the considered biomarkers (see.Appendix 1).
Several emulator architectures were tried on the training and validation data sets and the final choice was handpicked as a good tradeoff between high accuracy on the validation set (#1) and low computational runtime cost. We decided to utilize a twostage emulator architecture: First, the maximum conductances $\mathbf{x}$ – normalized to the range ${\mathbf{x}}_{i}\in [0.5,\text{}0.5]$ – are encoded using a first NN (${\mathrm{\Theta}}_{1}$) into a latent representation $\vartheta $. Second, this intermediate representation parameterizes a function ${f}_{\vartheta}:\mathbb{R}\to \mathbb{R}$ defined by a second NN (${\mathrm{\Theta}}_{2}$) that can be continuously evaluated to receive the emulated AP at time $t$. To help the second NN in computing the fast depolarization, a simple depolarization term $(\mathrm{tanh})$ is added to ${f}_{\vartheta}$. The parameters of this depolarization function are slope (d_{1}), offset (d_{2}), and amplitude (d_{3}), and are created by encoding the parameters through the first network, similar to the latent code. The AP approximated by the emulator is thus defined by
A schematic drawing of the emulator architecture is provided in 2. Splitting the network into two parts — one for encoding the parameters into a latent space and a second one for evaluating ${f}_{\vartheta}$ — allowed us to give the emulator enough complexity without markedly increasing the computational cost: in most cases, it is desirable to compute the whole AP in the entire range, for example $[0,\text{}1000]\text{}\mathrm{m}\mathrm{s}$, and not only at a single time step. Creating a single network that computes the mapping from maximum conductances to the transmembrane potential at a single time step (compare Figure 3) is either orders of magnitudes slower than encoding the parameters into a latent vector (only done once per AP) or would require to reduce the complexity of the network, which led to inaccurate emulations in the validation. The additional depolarization term was introduced to address the difficulty of fitting the depolarization phase during training and decreased the required training time substantially. Note that the mapping from maximum conductances $\mathbf{x}$ to depolarization parameters $\{{d}_{1},{d}_{2},{d}_{3}\}$ is also learned through ${\mathrm{\Theta}}_{1}$.
The exact architecture employed – chosen by a crossvalidation approach (see Validation) – comprised a first network (${\mathrm{\Theta}}_{1}$) of four fully connected layers of 256 neurons, each to encode the parameters into the latent vector $\vartheta \in {\mathbb{R}}^{256}$. This first network additionally generates the parameterization for the depolarization model $\{{d}_{1},{d}_{2},{d}_{3}\}$. The second network, computing the APs from the latent representation (${\mathrm{\Theta}}_{2}$), consisted of four fully connected layers of 64 neurons each. Exponential linear unit (ELU) activation functions Clevert et al., 2016 were used for all layers, except for the final nonlinear layer, which was modeled using a $\mathrm{tanh}$ activation function followed by a ${[1,1]}^{64}\to \mathbb{R}$ linear layer.
Training
Request a detailed protocolAlthough different regularization schemes such as variational losses (e.g. $\frac{\mathrm{d}{V}_{\mathrm{m}}}{\mathrm{d}t}$) were tried, the wealth of training data allowed us to define the training loss purely in terms of meansquarederror (MSE)
where ${f}_{{\vartheta}_{{\mathrm{\Theta}}_{1}}(\mathbf{x}),{\mathrm{\Theta}}_{2}}$ describes the output of the emulator using the current NN weights ${\mathrm{\Theta}}_{1}/{\mathrm{\Theta}}_{2}$, $\overline{X}$ refers to the current training batch and $X$ is the training data set containing both target APs ${\widehat{V}}_{\mathrm{m}}$ and corresponding maximum conductances $\mathbf{x}$. For the training, increasing batch sizes ($\overline{X}1250$ 1250 to 1800) were used, both in terms of the entire AP and considering only subsets over time $T\subset \mathbb{T}$ from $\frac{1}{16}\mathbb{T}$ to $\mathbb{T}$, resulting in a reduction of training time needed. The neural network was trained for a total number of 5000 epochs using the firstorder gradientbased algorithm ADAM Kingma and Ba, 2017. The training time was approximately $4\text{}\mathrm{h}$ on the GPU specified in Computational performance.
Validation
Request a detailed protocolA crossvalidation approach was used to quantify and compare the performance of various emulator architectures. The validation was based on 20% of the initial training data set (7976 pairs of maximum conductances and corresponding APs, see Training). For each pair, the emulated AP ${V}_{\mathrm{m}}$ was compared against the simulated AP ${\widehat{V}}_{\mathrm{m}}$, given the same maximum conductances. The mismatch was quantified by the rootmeansquared error (RMSE) defined as
The mismatch was also quantified in terms of AP biomarkers $\mathbf{b}\in {\mathbb{R}}^{N}$ (see AP biomarkers and abnormalities) and normalized maximum conductances $\mathbf{x}$ (see Time series fitting and estimation of maximum conductances and pharmacological parameters) in which case the RMSE was defined as
for $N$ samples.
Time series fitting and estimation of maximum conductances and pharmacological parameters
Request a detailed protocolTime series fitting is the basis for solving the inverse problem. To fit a given AP ${\widehat{V}}_{\mathrm{m}}$, defined on a subset of the trained domain $\hat{\mathbb{T}}\subset \mathbb{T}$, the first step was to choose a trial set of maximum conductances ${\mathbf{x}}_{0}$. Then, for the given trial set, the corresponding AP was emulated and the trial set of maximum conductances was iteratively updated by solving the following minimization problem:
where ${V}_{\mathrm{m}}(\mathbf{x},t):={f}_{{\theta}_{{\mathrm{\Theta}}_{1}}(\mathbf{x}),{\mathrm{\Theta}}_{2}}(t)$ is a shorthand for the emulator approximation function and t_{0} is a temporal offset parameter helping in fitting the exact depolarization timing. Here, ${\delta}_{[0.5,0.5]}(\mathbf{x})$ is the elementwise indicator function on the normalized feasible parameter space $[0.5,0.5]$. The minimization was done using ADAM (Kingma and Ba, 2017) combined with a projection on the feasible space.
To estimate maximum conductances for a given control AP, the control AP was fitted using the original maximum conductance values as initial guesses and priors: ${\mathbf{x}}_{0}=\widehat{\mathbf{x}}=0$. To estimate maximum conductances for a given drugged AP, the drugged AP was fitted using the maximum conductances estimated for the respective control AP as initial guesses and priors: ${\mathbf{x}}_{0}=\widehat{\mathbf{x}}={\mathbf{x}}^{c}$. The pharmacological parameters (scaling factors of control maximum conductances) were computed as elementwise ratios between the drugged and control maximum conductances $(\mathbf{s}{)}_{i}:=\frac{({\mathbf{G}}^{d}{)}_{i}}{({\mathbf{G}}^{c}{)}_{i}}$ but here, ${\mathbf{G}}^{d}$ and ${\mathbf{G}}^{c}$ are the nonnormalized maximum conductances, where $({\mathbf{G}}^{c}{)}_{i}>0$.
Since multiple sets of maximum conductances produced similarly good fits of the given AP, the parameter ${\lambda}_{{\mathbf{x}}_{0}}$ was introduced which minimizes the difference between original and control maximum conductances and between control and drugged maximum conductances, respectively. The value was chosen to be ${\lambda}_{{\mathbf{x}}_{0}}=10$ with respect to the highest accuracy found for the synthetic data set #2 that was generated for the evaluation (see Synthetic test data (#2/#3)).
Evaluation
The evaluation was performed for forward and inverse problems in pharmacological studies on synthetic and experimental data.
The raw experimental data were obtained without filtering but some filtering was applied before interfacing the data with the emulator. The APs contained a stimulus artifact between 0 and $1.5\text{}\mathrm{m}\mathrm{s}$ that was filtered as follows. For each of the last 10 consecutive APs, the transmembrane potential ${V}_{\mathrm{m}}$ closest to the end of the recorded time series was defined as the resting transmembrane potential $RMP$ and ${V}_{\mathrm{m}}(t)=RMP$ was set for $t\in [0,\text{}1.75]\mathrm{m}\mathrm{s}$. Then, the APs were resampled at $100\text{}\mathrm{k}\mathrm{H}\mathrm{z}$ and a lowpass filtering was performed with a secondorder butterworth filter (cutoff at $2.5\text{}\mathrm{k}\mathrm{H}\mathrm{z}$) to reduce the highfrequency noise of the signal. Finally, the filtered APs were averaged and the averaged AP was again resampled at $1\text{}\mathrm{k}\mathrm{H}\mathrm{z}$ for $t\in [15,\text{}1000]\mathrm{m}\mathrm{s}$ (repolarization) and $100\text{}\mathrm{k}\mathrm{H}\mathrm{z}$ in $t\in [0,\text{}15)\mathrm{m}\mathrm{s}$ (depolarization). An example comparison of a raw and a filtered averaged AP is given in Figure 4.
Computational performance
Request a detailed protocolThe simulation of a single AP (see Simulator) sampled at a resolution of 20 kH_{Z} took 293 s on one core of a AMD Ryzen Threadripper 2990 WX (clock rate: 3.0 GH_{Z}) in CARPentry. Adaptive timestep solver of variable order, such as implemented in Myokit (Clerx et al., 2016), can significantly lower the simulation time (30 s for our setup) by using small step sizes close to the depolarization (phase 0) and increasing the time step in all other phases. The emulation of a steady state AP sampled at a resolution of 20 KH_{Z} for $t\in [10,1000]\text{}\mathrm{m}\mathrm{s}$ took 18.7 ms on a AMD Ryzen 7 3800 X (clock rate: 3.9 GH_{Z}) and 1.2 ms on a Nvidia A100 (Nvidia Corporation, USA), including synchronization and data copy overhead between CPU and GPU.
The amount of required beats to reach the steady state of the cell in the simulator has a major impact on the runtime and is not known apriori. On the other hand, both simulator and emulator runtime linearly depends on the time resolution, but since the output of the emulator is learned, the time resolution can be chosen at arbitrarily without affecting the AP at the sampled times. This makes direct performance comparisons between the two methodologies difficult. To still be able to quantify the speedup, we ran Myokit using 100 beats to reach steady state, taking 3.2 s of simulation time. In this scenario, we witnessed a speedup of 171 and $2\cdot {10}^{3}$ of our emulator on CPU and GPU, respectively (again including synchronization and data copy overhead between CPU and GPU in the latter case). Note that both methods are similarly expected to have a linear parallelization speedup across multiple cells.
For the inverse problem, we parallelized the problem for multiple cells and keep the problem on the GPU to minimize the overhead, achieving emulations (including backpropagation) that run in 120 s per AP at an average temporal resolution of 2 KH_{Z}. We consider this the peak performance which will be necessary for the inverse problem in Inverse problem based on synthetic data.
Forward problem
Request a detailed protocolThe emulator evaluation for the forward problem, i.e. to find the drugged AP for given pharmacological parameters, was only performed on synthetic data since maximum conductances were not available experimentally. The maximum conductances of data sets #2 and #3 were used to consider data with normal APs and with abnormal APs exhibiting EADs. Pharmacological parameters are not inputs of the emulator but drugged maximum conductances that were computed as control maximum conductances scaled by the given pharmacological parameters (see Synthetic test data (#2/#3)). These were used to emulate drugged APs. The RMSE was used to quantify the mismatch between the emulated and the ground truth AP.
Inverse problem
Request a detailed protocolThe emulator evaluation for the inverse problem, that is to find the pharmacological parameters for given control and drugged APs (through optimization), was performed on both synthetic and experimental data. When using synthetic data, the data set #2 was used including data with normal APs. First, control and drugged maximum conductances were estimated based on control and drugged APs and then, pharmacological parameters were computed as ratios of drugged and control conductances (see Time series fitting and estimation of maximum conductances and pharmacological parameters). The mismatch between estimated and ground truth maximum conductances were quantified using the error that is defined through
where $\mathbf{x}$ and $\widehat{\mathbf{x}}$ are the normalized estimated and ground truth maximum conductances (see Architecture). Similarly, we computed the mismatch between estimated and ground truth scaling factor vectors ($\mathbf{S}$ and $\widehat{\mathbf{S}}$ respectively) as
When using experimental data, maximum conductances and pharmacological parameters were estimated in the same way but due to a lack of experimental maximum conductances, the mismatch between estimated and ground truth values could not be quantified. Instead, the estimated pharmacological parameters were compared with distributions computed from data published within the CiPA initiative (Li et al., 2017; Chang et al., 2017) (CiPA distributions). The data set (https://github.com/FDA/CiPA/tree/master/Hill_fitting/results; Chang and Li, 2017) includes 2000 IC_{50} values and Hill coefficients for each drug and for up to seven targets (${I}_{Na}$, ${I}_{NaL}$, ${I}_{CaL}$, ${I}_{to}$, ${I}_{Kr}$, ${I}_{Ks}$, ${I}_{K1}$). The poreblock model (Brennan et al., 2009) was used to obtain the corresponding scaling factors.
Results
Evaluation
Forward problem
The emulator evaluation for the forward problem was only done on synthetic data and both data sets #2 and #3 (see Synthetic test data (#2/#3)) were used to analyze the solution accuracy for normal and abnormal APs exhibiting EADs.
The data set #2 was used first and Figure 5 illustrates the distribution of RMSEs between emulated and ground truth drugged APs. In total 10^{4} APs were emulated in 0.6 s. The average RMSE over all APs was 0.47 mV and only for a few APs the RMSE was >1 mV with 1.5 mV being the maximum. Largest mismatches were located in the phases 0 and 3 of the AP. While the mismatches in phase 3 were simply a result of imperfect emulation, the mismatches in phase 0 were a result of the difficulty in matching the depolarization time exactly.
Figure 6 shows the distribution of biomarker mismatches between emulated and ground truth drugged APs. The low RMSEs between the APs translated into low RMSEs between the AP biomarkers. Likewise, the difficulty in exactly matching the depolarization time leads to elevated errors and more outliers in the biomarkers influenced by the depolarization phase ($TP$ and $dVmMax$).
The data set #3 was used second and Appendix 3 shows all emulated APs, both containing the EAD and nonEAD cases. The emulation of all 950 APs took 0.76 s on the GPU specified in Training We show the emulation of all maximum conductances and the classification of the emulation. The comparison with the actual EAD classification (based on the criterion outlined in Appendix 1) results in truepositive (EAD both in the simulation and emulation), falsenegative (EAD in the simulation, but not in the emulation), falsepositive (EAD in the emulation, but not in the simulation) and truenegative (no EAD both in the emulation and simulation). The emulations achieved 72.5% sensitivity (EAD cases correctly classified) and 94.9% specificity (nonEAD cases correctly classified), with an overall acurracy of 90.8% (total samples correctly classified). A substantial amount of wrongly classified APs showcase a notable proximity to the threshold of manifesting EADs. Figure 7 illustrates the distribution of RMSEs in the EAD APs between emulated and ground truth drugged APs. The average RMSE over all EAD APs was 14.5 mV with 37.1 mV being the maximum. Largest mismatches were located in phase 3 of the AP, in particular in emulated APs that did not fully repolarize.
Inverse problem based on synthetic data
The emulator evaluation for the inverse problem was first done using synthetic data (data set #2, see Synthetic test data (#2/#3)). For this, we minimized (3) by using ADAM (Kingma and Ba, 2017) with no batching for 10^{4} iterations for all 100 cardiomyocytes times 100 drugs (i.e. 10^{4} APs), resulting in a total of 10^{8} emulations, taking approximately 3.5 h on the GPU specified in Training. Control and drugged APs could be fitted with an average RMSE of 0.8 mV. Largest mismatches were located in phase 0 and 3 of the AP for the reasons given above (see Forward problem). Figure 8 shows the distribution of the errors between the estimated and the ground truth maximum conductances and pharmacological parameters. For both the maximum conductances (RMSE $\le 0.18$) and the related pharmacological parameters (RMSE $\le \phantom{\rule{thinmathspace}{0ex}}0.22$), the errors were closely distributed around zero. However, the RMSEs increased from the control maximum conductance over the drugged maximum conductance to the pharmacological parameters and there were distinctive differences among maximum conductances and related pharmacological parameters with the smallest for G_{Kr} and the largest for G_{Ks}.
Inverse problem based on experimental data
The emulator evaluation for the inverse problem was then done on experimental data (Experimental data (#4)). Similar to the synthetic inverse problem, we optimized (3), this time for $5\cdot {10}^{4}$ and $2.5\cdot {10}^{4}$ epochs for control and drug APs for each drug sample. The total computational time spent on this task was 2.8 hr on the GPU specified in Training.
Figure 9 shows the fitted and the ground truth APs for all drugs. Control and drugged APs could be fitted with average RMSEs shown in Table 4. The largest mismatch was located in phase 0 for most APs, which was the result from an imperfect matching of the exact depolarization timing.
Figure 10 compares the estimated pharmacological parameters and the CiPA distributions. Pharmacological parameters that fell into the range spanned by $\mu \pm \left(0.15+\sigma \right)$ of the CiPA distribution, where $\mu ,\sigma $ are the distribution’s mean and standard deviation respectively, were classified as successfully estimated while the others were classified as unsuccessfully estimated. In total, 50% of the pharmacological parameters could be estimated succesfully and while all pharmacological parameters related to G_{Ks} could be successfully estimated, unsuccessfully estimated parameters were found across all maximum conductances, in particular related to G_{K1} for which all pharmacological parameters could not be successfully estimated (Table 5).
Discussion
NN emulation of the human ventricular cardiomyocyte AP was introduced and the applicability in pharmacological studies was investigated.
Evaluation
The evaluated NN emulator showed highly accurate AP emulations for the forward problem on synthetic data. High accuracy was found in normal APs in data set #2 (average RMSE was 0.47 mV; Figure 5) and to a lesser extent also in abnormal APs exhibiting EADs: of the emulated EAD APs, 72.5% exhibited alignment with the abnormality, and the substantial majority of the remaining APs demonstrated pronounced proximity, while the average RMSE among the EAD APs was 14.5 mV (Figure 7). In comparison, the normal APs exhibiting no EAD in data set #3 could be reconstructed with an average RMSE of 4.02 mV and the overall accuracy of classifying EAD on emulated APs was 90.8%. Increasing the amount of training data within the relevant range could lead to further enhancements in accuracy for abnormal APs. Nevertheless, this observation demonstrates that the emulator is also capable of accounting for discontinuities of the response surface. This is particularly useful in pharmacological studies and a key advantage over existing emulation approaches (Chang et al., 2015; Johnstone et al., 2016; Coveney and Clayton, 2018; Ghosh, 2018; Rasmussen, 2019; Coveney et al., 2021).
The emulator was further evaluated for the inverse problem on synthetic and also on experimental data. Maximum conductances and related pharmacological parameters could be widely estimated with high accuracy on synthetic data (RMSE $\le 0.18$ and $\le 0.21$ for all maximum conductances and pharmacological parameters, respectively; Figure 8).
The RMSEs increased from the control maximum conductance over the drugged maximum conductance to the pharmacological parameters which may be because the estimation of drugged maximum conductances depends on the control maximum conductances and the computation of pharmacological parameters depends on both control and drugged maximum conductances (see Time series fitting and estimation of maximum conductances and pharmacological parameters) allowing that errors can propagate and amplify. Distinctive differences were observed among the maximum conductances and related pharmacological parameters and the largest RMSEs were found for G_{Ks} throughout. This can be attributed to various degrees of parameter identifiability Sarkar and Sobie, 2010; Zaniboni et al., 2010; Groenendaal et al., 2015; Jaeger et al., 2019; Whittaker et al., 2020 and the results agree with the GSA that indicates almost nonidentifiability of G_{Ks} (Appendix 2).
Larger inaccuracies were found in the inverse problem solutions on experimental data (Figure 10, Table 4). The first reason may be low parameter identifiability and we want to highlight inaccuracies in estimating the pharmacological parameters related to G_{Kr}, P_{Ca}, and G_{NaL} when the hERG channel was blocked in parallel to the Cav1.2 channel (verapamil) or in parallel to both the Cav1.2 and the Nav1.5late channel (terfenadine). The hERG channel block (prolongation of the AP), and the Cav1.2 and Nav1.5late channel block (shortening of the AP) are known to have opposite effects on the AP Orvos, 2019. At the given drug concentrations, these effects were apparently counterbalancing, which resulted in negligible changes of the AP (Figure 9). This situation made the estimation of pharmacological parameters very challenging and led to particularly large inaccuracies for terfenadine. The accurate estimations of the pharmacological parameters related to GKs are surprising at first in light of the almost nonidentifiability. This was due to the combination of two factors: (1) different from synthetic data, the drugs at the given concentrations did not affect the corresponding KCNQ1MinK channel and (2) the difference between control and drugged maximum conductance was weakly enforced to be minimal (see Time series fitting and estimation of maximum conductances and pharmacological parameters) which leads to almost no difference in nonidentifiable parameters and hence, to a pharmacological parameter of one. The accuracy will likely be much lower in drugs that affect the KCNQ1MinK channel.
The second and probably main reason for the inaccuracies may be the fact that the data were collected in small tissue preparations, whereas the emulator was trained on data generated by a simulator that represents single cardiomyocytes. APs in small tissue preparations are slightly different from those in single cardiomyocytes. Differences can arise from electrotonic coupling and the mixture of cells including fibroblasts that are able to modify the EP (Kohl and Gourdie, 2014; Mayer et al., 2017; Hall et al., 2021). This can hamper the fitting of the APs and consequently, the estimation of the maximum conductances and pharmacological parameters.
Emulator
The presented NN emulator enables a massive speedup compared to regular simulations and the evaluation for the forward problem on synthetic data showed also highly accurate AP emulations. Cardiomyocyte EP models are already very quick to evaluate in the scale of seconds (see Computational performance), but the achieved runtime of emulations allows to solve time consuming simulation protocols markedly more efficient. One such scenario is the presented inverse maximum conductance estimation problem (see Inverse problem based on synthetic data and Inverse problem based on experimental data), where for estimating maximum conductances of a single AP, we need to emulate the steady state AP at least several hundred times as part of an optimization procedure. Further applications include the probabilistic use of cardiomyocyte EP models with uncertainty quantification (Chang et al., 2017; Johnstone et al., 2016) where thousands of samples of parameters are potentially necessary to compute a distribution of the steadystate properties of subsequent APs, and the creation of cell populations (Muszkiewicz et al., 2016; Gemmell et al., 2016; Britton et al., 2013). In addition to the aforementioned strengths of the presented NN emulator, some further valuable features are worth mentioning that arise from the continuous nature of the emulator. First, the AP can be emulated and fitted at any desired resolution. Second, timing offsets, for example between stimuli in the data to be fitted and the training data, can be accounted for using t_{0} in (3) without retraining. Last but not least, the transmembrane potential gradient $\frac{\mathrm{d}{V}_{\mathrm{m}}}{\mathrm{d}t}$, recently highlighted in terms of proarrhythmic potential prediction (Jeong et al., 2022), can be continuously derived and is not dependent on the temporal discretization.
Limitations and future work
Some limitations have to be considered. First, the emulator has only maximum conductances as inputs. Although these explain much of the AP variability seen between cardiomyocytes (Britton et al., 2013; Muszkiewicz et al., 2016), the inclusion of parameters related to the channel kinetics might enable a more detailed consideration of drug effects in pharmacological studies. In general, the question of complete parameter identifiability utilizing only APs remains an open challenge (Zaniboni et al., 2010). Channel kinetics determine the contribution of the corresponding current to the AP generation in different phases and can thus also modulate drug effects, but were however neglected, as the expansion of the input space might be unsuitable for solving the inverse problem when only AP data are used. Second, the interaction between drugs and their targets is solely captured through scaling of the related maximum conductance at control, which is mostly adequate but in fact an oversimplification (Brennan et al., 2009). The interaction can be dependent on time, voltage, and channel state, which requires the use of Markov models with many more parameters (Brennan et al., 2009; Li et al., 2017; Lei et al., 2024). Again, this expands the input space and might be unsuitable for solving the inverse problem when only AP data are used. Moreover, drugs that are applied over a longer period of time can also cause modifications of the maximum conductances through changes in gene expression (Shim et al., 2023). This requires attention to avoid misinterpretations of found blocking or enhancement effects, for example by estimating the control maximum conductances again after a washout procedure. Third, the inverse problem was only solved for AP data obtained from one single stimulation protocol. Johnstone et al., 2016 have shown that the usage of AP data obtained from various stimulation protocols can improve the parameter identifiability and thus, the accuracy of parameter estimates. To be able to use those data in the presented approach, the pacing cycle length must be included as additional input in the emulator and the emulator may be trained on more than the last AP of the pacing series. This would also allow to capture alternans. Last but not least, the number of drugs and concentrations considered in the inverse problem on experimental data poses a limitation. The ultimate goal is to have a tool that provides highly accurate solutions for drugs with different targets and concentrations. To this end, analyses must be extended by data obtained from a series of available and well characterized drugs. The data should be collected in single cardiomyocytes in order to minimize the discussed inaccuracies that stem from the use of tissue preparation data. Further simulation studies of the inverse problem on tissue slabs versus cardiomyocyte EP model could be integrated to assess the impact of differences in setups. This should be addressed in future work. Of note, the presented approach can also be straightforwardly applied to other transients, for example intracellular [Ca^{2+}] or sarcomere length.
Conclusion
This paper introduced NN emulation of the human ventricular cardiomyocyte AP and tested its applicability in pharmacological studies. The computational cost of the NN emulator was compared to that of the simulator, revealing a massive speedup of more than 1e^{3}. The accuracy of solving the forward problem on synthetic data was found to be high for normal APs and this hold mostly true for abnormal APs exhibiting EADs. This advantage distinguishes our novel approach from existing emulation methods. While larger inaccuracies were observed when utilizing experimental data – a limitation thoroughly discussed and particularly inherent to the fact that small tissue preparations were studied while the emulator was trained on single cardiomyocyte data – the accuracy of solving the inverse problem on synthetic data remained high. Collectively, these findings underscore the potential of NN emulators in improving the efficiency of future quantitative systems pharmacology studies.
Appendix 1
AP biomarkers and abnormalities
The AP biomarkers used in this study were selected such that the key characteristics of the depolarization and the repolarization phase can be quantified. They include $RMP$ (resting transmembrane potential measured just before stimulation), $dVmMax$ (maximum transmembrane potential slope during the upstroke), $Peak$ (peak transmembrane potential at the end of the upstroke), $TP$ (timetopeak from the stimulus until $Peak$ is reached), $AP{D}_{x}$ (AP duration at $x\in \{30,40,50,60,70,80,90\}$% repolarization relative to the AP amplitude ($PeakRMP$) measured from the instant of $dVmMax$), and $Tr{i}_{9040}$ (triangulation defined as the difference between $AP{D}_{90}$ and $AP{D}_{40}$ Britton et al., 2017).
For the GSA (see.3) and the creation of synthetic data (see Synthetic test data (#2/#3)), AP abnormalities were also considered which included depolarization abnormalities, repolarization abnormalities, and alternans. Depolarization abnormalities were defined as an upstroke peak below 0 mV and an AP that does not reach 0 mV before 100 ms after stimulation (Passini et al., 2017). Repolarization abnormalities were defined as a transmembrane potential rate of rise of more than $0.01\text{}\frac{\mathrm{mV}}{\mathrm{ms}}$ from 150 ms after the upstroke peak onwards (representative of early afterdepolarizations) and as a transmembrane potential that does not fall below 40 mV (Passini et al., 2017) (representative of repolarization failure). Alternans were defined as $AP{D}_{90}$ difference of more than 5 ms between two consecutive APs (Morotti et al., 2021).
Appendix 2
Global sensitivity analysis
A variancebased Sobol’ global sensitivity analysis (GSA) (Sobol′, 2001) was performed on the simulator to quantify the sensitivities of the maximum channel conductances (inputs) with respect to the AP biomarkers (outputs; see AP biomarkers and abnormalities). This informed the decision on which inputs to consider in the emulator. Furthermore, it was used for the interpretation of the solutions of the inverse problem since parameters that are insensitive with respect to the outputs indicate nonidentifiability (Guillaume et al., 2019).
The maximum conductances used for building the model population in Tomek et al., 2019 were considered and Saltelli’s sampling scheme (Saltelli, 2002) was applied with $N=1024$ to generate 20,480 input samples with values between 50% and 150% of the original values. Simulations were performed for each input sample and biomarker values (see AP biomarkers and abnormalities) derived from the last AP were used for the analysis. However, data were excluded if not all biomarkers could be determined or abnormalities (see AP biomarkers and abnormalities) were detected in the last two consecutive APs. Firstorder (S1) and totaleffect (ST) Sobol’ sensitivity indices were computed using the Saltelli method (Homma and Saltelli, 1996; Saltelli, 2002). This requires outputs for each input sample and to take this into account, excluded outputs were assigned the mean values of included outputs. The SALibSensitivity Analysis Library (Herman and Usher, 2017) was used for the GSA.
The GSA could be performed on the data of all input samples as no data were excluded. S1 and ST were mostly very similar which indicates only little interactions among the maximum conductivities relative to the AP biomarkers (9). The only exception was $TP$. As was to be expected, the analysis underlines the predominant relative sensitivity of G_{Na} with respect to biomarkers of the depolarization phase, the predominant relative sensitivity of G_{NaL}, G_{Kr}, and G_{NCX} with respect to biomarkers of the repolarization phase and the predominant relative sensitivity of G_{K1}, P_{NaK} to the resting transmembrane potential. However, G_{Ks} has a negligible relative sensitivity to all biomarkers. This indicates almost nonidentifiability.
Appendix 3
EAD classification
Data availability
The trained emulator is available as a python package from https://github.com/thomgrand/cardiomyocyte_emulator (copy archived at Thomas, 2024; the code is licensed under AGPLv3, see https://www.gnu.org/licenses/agpl3.0.en.html for details). The trained emulator is provided as a Python package, heavily utilizing PyTorch (Paszke, 2019) for the neural network execution, allowing it to be executed on both CPUs and NVidia GPUs. Additionally, NumPy (Harris et al., 2020) interfaces are provided for easy interfacing with other libraries. Note that the provided repository is a reimplementation of the codebase used in this study and thus may deviate in performance or runtime. The training data set (data set #1) is available on Zenodo using the link https://zenodo.org/records/10640339.
References

Evolution of mathematical models of cardiomyocyte electrophysiologyMathematical Biosciences 334:108567.https://doi.org/10.1016/j.mbs.2021.108567

Preclinical cardiac safety assessment of pharmaceutical compounds using an integrated systemsbased computer model of the heartProgress in Biophysics and Molecular Biology 90:414–443.https://doi.org/10.1016/j.pbiomolbio.2005.06.006

Multiscale modelling of druginduced effects on cardiac electrophysiological activityEuropean Journal of Pharmaceutical Sciences 36:62–77.https://doi.org/10.1016/j.ejps.2008.09.013

Myokit: a simple interface to cardiac cellular electrophysiologyProgress in Biophysics and Molecular Biology 120:100–114.https://doi.org/10.1016/j.pbiomolbio.2015.12.008

The comprehensive in vitro proarrhythmia assay (CiPA) initiative  update on progressJournal of Pharmacological and Toxicological Methods 81:15–20.https://doi.org/10.1016/j.vascn.2016.06.002

Using cardiac ionic cell models to interpret clinical dataWIREs Mechanisms of Disease 13:e1508.https://doi.org/10.1002/wsbm.1508

Fitting two human atrial cell models to experimental data using Bayesian history matchingProgress in Biophysics and Molecular Biology 139:43–58.https://doi.org/10.1016/j.pbiomolbio.2018.08.001

Bayesian calibration of electrophysiology models using restitution curve emulatorsFrontiers in Physiology 12:693015.https://doi.org/10.3389/fphys.2021.693015

An evaluation of 30 clinical drugs against the comprehensive in vitro proarrhythmia assay (CiPA) proposed ion channel panelJournal of Pharmacological and Toxicological Methods 81:251–262.https://doi.org/10.1016/j.vascn.2016.03.009

OpenCOR: a modular and interoperable approach to computational biologyFrontiers in Physiology 6:26.https://doi.org/10.3389/fphys.2015.00026

Rabbitspecific computational modelling of ventricular cell electrophysiology: Using populations of models to explore variability in the response to ischemiaProgress in Biophysics and Molecular Biology 121:169–184.https://doi.org/10.1016/j.pbiomolbio.2016.06.003

A novel computational model of the human ventricular action potential and Ca transientJournal of Molecular and Cellular Cardiology 48:112–121.https://doi.org/10.1016/j.yjmcc.2009.09.019

Cellspecific cardiac electrophysiology modelsPLOS Computational Biology 11:e1004242.https://doi.org/10.1371/journal.pcbi.1004242

Introductory overview of identifiability analysis: A guide to evaluating whether you have the right type of data for your modeling purposeEnvironmental Modelling & Software 119:418–432.https://doi.org/10.1016/j.envsoft.2019.07.007

Complex relationship between cardiac fibroblasts and cardiomyocytes in health and diseaseJournal of the American Heart Association 10:e019338.https://doi.org/10.1161/JAHA.120.019338

SALib: an opensource python library for sensitivity analysisThe Journal of Open Source Software 2:97.https://doi.org/10.21105/joss.00097

Importance measures in global sensitivity analysis of nonlinear modelsReliability Engineering & System Safety 52:1–17.https://doi.org/10.1016/09518320(96)000026

Adaptive activation functions accelerate convergence in deep and physicsinformed neural networksJournal of Computational Physics 404:109136.https://doi.org/10.1016/j.jcp.2019.109136

Uncertainty and variability in models of the cardiac action potential: Can we build trustworthy models?Journal of Molecular and Cellular Cardiology 96:49–62.https://doi.org/10.1016/j.yjmcc.2015.11.018

Fibroblastmyocyte electrotonic coupling: does it occur in native cardiac tissue?Journal of Molecular and Cellular Cardiology 70:37–46.https://doi.org/10.1016/j.yjmcc.2013.12.024

Improving cardiomyocyte model fidelity and utility via dynamic electrophysiology protocols and optimization algorithmsThe Journal of Physiology 594:2525–2536.https://doi.org/10.1113/JP270618

Neural network differential equations for ion channel modellingFrontiers in Physiology 12:708944.https://doi.org/10.3389/fphys.2021.708944

The impact of uncertainty in hERG binding mechanism on in silico predictions of druginduced proarrhythmic riskBritish Journal of Pharmacology 181:987–1004.https://doi.org/10.1111/bph.16250

Improving the in silico assessment of proarrhythmia risk by combining hERG (human EtheràgogoRelated Gene) channeldrug binding kinetics and multichannel pharmacologyCirculation. Arrhythmia and Electrophysiology 10:e004628.https://doi.org/10.1161/CIRCEP.116.004628

Assessment of an in silico mechanistic model for proarrhythmia risk prediction under the CiPA InitiativeClinical Pharmacology and Therapeutics 105:466–475.https://doi.org/10.1002/cpt.1184

Complex restitution behavior and reentry in a cardiac tissue model for neonatal micePhysiological Reports 5:e13449.https://doi.org/10.14814/phy2.13449

Prediction of Thorough QT study results using action potential simulations based on ion channel screensJournal of Pharmacological and Toxicological Methods 70:246–254.https://doi.org/10.1016/j.vascn.2014.07.002

Uncertainty and variability in computational and mathematical models of cardiac physiologyThe Journal of Physiology 594:6833–6847.https://doi.org/10.1113/JP271671

Emulating dynamic nonlinear simulators using Gaussian processesComputational Statistics & Data Analysis 139:178–196.https://doi.org/10.1016/j.csda.2019.05.006

Variability in cardiac electrophysiology: Using experimentallycalibrated populations of models to move beyond the single virtual physiological human paradigmProgress in Biophysics and Molecular Biology 120:115–127.https://doi.org/10.1016/j.pbiomolbio.2015.12.002

Probabilistic sensitivity analysis of complex models: a bayesian approachJournal of the Royal Statistical Society Series B 66:751–769.https://doi.org/10.1111/j.14679868.2004.05304.x

The virtual assay software for human in silico drug trials to augment drug cardiac testingJournal of Computational Science 52:101202.https://doi.org/10.1016/j.jocs.2020.101202

ConferencePyTorch: An Imperative Style, HighPerformance Deep Learning LibraryAdvances in Neural Information Processing Systems.

The openCARP simulation environment for cardiac electrophysiologyComputer Methods and Programs in Biomedicine 208:106223.https://doi.org/10.1016/j.cmpb.2021.106223

A practical algorithm for solving dynamic membrane equationsIEEE Transactions on BioMedical Engineering 25:389–392.https://doi.org/10.1109/TBME.1978.326270

Making best use of model evaluations to compute sensitivity indicesComputer Physics Communications 145:280–297.https://doi.org/10.1016/S00104655(02)002801

Regression analysis for constraining free parameters in electrophysiological models of cardiac cellsPLOS Computational Biology 6:e1000914.https://doi.org/10.1371/journal.pcbi.1000914

A quantitative systems pharmacology perspective on the importance of parameter identifiabilityBulletin of Mathematical Biology 84:39.https://doi.org/10.1007/s11538021009825

Predicting individualspecific cardiotoxicity responses induced by tyrosine kinase inhibitorsFrontiers in Pharmacology 14:1158222.https://doi.org/10.3389/fphar.2023.1158222

Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimatesMathematics and Computers in Simulation 55:271–280.https://doi.org/10.1016/S03784754(00)002706

Translational models and tools to reduce clinical trials and improve regulatory decision making for QTc and proarrhythmia risk (ICH E14/S7B updates)Clinical Pharmacology and Therapeutics 109:319–333.https://doi.org/10.1002/cpt.2137

Atrial cell action potential parameter fitting using genetic algorithmsMedical & Biological Engineering & Computing 43:561–571.https://doi.org/10.1007/BF02351029

SoftwareCardiomyocyte_Emulator, version swh:1:rev:26ba43c619c0acf2fc853326aa4cd0578b786d6aSoftware Heritage.

Mechanistic modelinformed proarrhythmic risk assessment of drugs: review of the “CiPA” initiative and design of a prospective clinical validation studyClinical Pharmacology and Therapeutics 103:54–66.https://doi.org/10.1002/cpt.896

Solvers for the cardiac bidomain equationsProgress in Biophysics and Molecular Biology 96:3–18.https://doi.org/10.1016/j.pbiomolbio.2007.07.012

Early afterdepolarizations and cardiac arrhythmiasHeart Rhythm 7:1891–1899.https://doi.org/10.1016/j.hrthm.2010.09.017

Calibration of ionic and cellular cardiac electrophysiology modelsWiley Interdisciplinary Reviews. Systems Biology and Medicine 12:e1482.https://doi.org/10.1002/wsbm.1482
Article and author information
Author details
Funding
Austrian Science Fund (ERANET I 4652B)
 Christoph M Augustin
National Institutes of Health (R01HL158667)
 Christoph M Augustin
German Research Foundation (Walter Benjamin Fellowship No. 468256475)
 Alexander Jung
HUNREN Hungarian Research Network
 Norbert Jost
 András Varró
Wellcome Trust (Senior Research Fellowship No. 212203/Z/18/Z)
 Gary R Mirams
National Research Development and Innovation Office (No. 142738)
 Norbert Jost
 András Varró
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication. For the purpose of Open Access, the authors have applied a CC BY public copyright license to any Author Accepted Manuscript version arising from this submission.
Acknowledgements
This research was funded by the Austrian Science Fund (FWF), the German Research Foundation (DFG), the National Institutes of Health, the National Research, Development and Innovation Office (NRDI), the Hungarian Research Network (HUNREN), and the Wellcome Trust. For the purpose of open access, the authors have applied a CCBY public copyright licence to any Author Accepted Manuscript version arising from this submission.
Ethics
Data set #4 contains data of small right ventricular trabeculae and papillary tissue preparations that were obtained from healthy human hearts (see Orvos et al., 2019). Hearts were obtained from organ donors whose hearts were explanted to obtain pulmonary and aortic valves for transplant surgery. Before cardiac explantation, organ donors did not receive medication apart from dobutamine, furosemide, and plasma expanders. Proper consent was obtained for use of each individual's tissue for experimentation. The investigations conform to the principles of the Declaration of Helsinki. Experimental protocols were approved by the University of Szeged and National Scientific and Research Ethical Review Boards (No. 5157/1997 OEj and 49910/20101018EKU [339/PI/010]).
Version history
 Preprint posted:
 Sent for peer review:
 Reviewed Preprint version 1:
 Reviewed Preprint version 2:
 Version of Record published:
Cite all versions
You can cite all versions using the DOI https://doi.org/10.7554/eLife.91911. This DOI represents all versions, and will always resolve to the latest one.
Copyright
© 2023, Grandits et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 397
 views

 55
 downloads

 1
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Cell Biology
 Developmental Biology
Chondrocyte columns, which are a hallmark of growth plate architecture, play a central role in bone elongation. Columns are formed by clonal expansion following rotation of the division plane, resulting in a stack of cells oriented parallel to the growth direction. In this work, we analyzed hundreds of Confetti multicolor clones in growth plates of mouse embryos using a pipeline comprising 3D imaging and algorithms for morphometric analysis. Surprisingly, analysis of the elevation angles between neighboring pairs of cells revealed that most cells did not display the typical stacking pattern associated with column formation, implying incomplete rotation of the division plane. Morphological analysis revealed that although embryonic clones were elongated, they formed clusters oriented perpendicular to the growth direction. Analysis of growth plates of postnatal mice revealed both complex columns, composed of ordered and disordered cell stacks, and small, disorganized clusters located in the outer edges. Finally, correlation between the temporal dynamics of the ratios between clusters and columns and between bone elongation and expansion suggests that clusters may promote expansion, whereas columns support elongation. Overall, our findings support the idea that modulations of division plane rotation of proliferating chondrocytes determines the formation of either clusters or columns, a multifunctional design that regulates morphogenesis throughout pre and postnatal bone growth. Broadly, this work provides a new understanding of the cellular mechanisms underlying growth plate activity and bone elongation during development.

 Cell Biology
 Immunology and Inflammation
NLRP3 is an inflammasome seeding pattern recognition receptor activated in response to multiple danger signals which perturb intracellular homeostasis. Electrostatic interactions between the NLRP3 polybasic (PB) region and negatively charged lipids on the transGolgi network (TGN) have been proposed to recruit NLRP3 to the TGN. In this study, we demonstrate that membrane association of NLRP3 is critically dependant on Sacylation of a highly conserved cysteine residue (Cys130), which traps NLRP3 in a dynamic Sacylation cycle at the Golgi, and a series of hydrophobic residues preceding Cys130 which act in conjunction with the PB region to facilitate Cys130 dependent Golgi enrichment. Due to segregation from Golgi localised thioesterase enzymes caused by a nigericin induced breakdown in Golgi organisation and function, NLRP3 becomes immobilised on the Golgi through reduced deacylation of its Cys130 lipid anchor, suggesting that disruptions in Golgi homeostasis are conveyed to NLRP3 through its acylation state. Thus, our work defines a nigericin sensitive Sacylation cycle that gates access of NLRP3 to the Golgi.