We model the electrophysiological behavior of cardiac tissue using the monodomain model (Clayton et al., 2011; Sahli Costabal et al., 2018a). This model's main variable is the transmembrane potential ϕ, the difference between the intra- and extra-cellular potentials. The transmembrane potential is governed by a reaction-diffusion equation (Krishnamoorthi et al., 2014),

Here, we introduce the source term f^ϕ which represents the ionic currents across the cell membrane and the conductivity tensor D that accounts for a fast signal propagation of D^|| parallel to the cardiac muscle fiber direction f and a slow signal propagation D^⊥ perpendicular to it (Clerc, 1976; Plank et al., 2008; Goktepe and Kuhl, 2009),

In general, the ionic currents f^ϕ are functions of the transmembrane potential ϕ and a set of state variables q(ϕ) (Wong et al., 2013; Lee et al., 2016), f^ϕ(ϕ, q(ϕ);t) where the state variables themselves are governed by ordinary differential equations, q˙=g(ϕ,q(ϕ);t). The number of currents and state variables determines the complexity of the cell model and varies for different cell types. For human Purkinje fiber cells, we adopt the Stewart model (Stewart et al., 2009) which tracks 14 ionic currents using 20 state variables

To represent electrophyiological behavior of human ventricular cardiomyocytes, we adopt the O'Hara-Rudy model (O'Hara et al., 2011) with a minor modification (Priest et al., 2016) of the fast sodium current I_Na (ten Tusscher et al., 2004). Studies have shown that this I_Na substitution leads to a physiological conduction velocity restitution behavior, with a minor impact on the action potential behavior (Elshrif and Cherry, 2014). The resulting model tracks 15 ionic currents defined through a total of 39 state variables,

To incorporate drug effects into our multiscale models, we selectively block the relevant ionic currents in the Purkinje and cardiomyocyte cell models (Sahli Costabal et al., 2018b). These blocks are informed by experimental patch-clamp experiments that study the fractional blockage β of different ion channels at varying drug concentrations (McMillan et al., 2017). We implement these fractional blockings using fitted Hill-type equations of the form,

which are characterized by the exponent h and the concentration IC50h required to achieve a 50% current block. To apply a specific drug, we select a desired concentration C, calculated the fractional blockage β_ion for each considered ion channel, and scale the corresponding ion current channels by the fractional blockage [1 − β],

Using the male and female multiscale cardiac electrophysiology models, we develop two sex-specific arrhythmogenic risk classifiers based on drug- and dose-specific ion channel blockage. Given the high computational cost of evaluating arrhythmogenesis for a single full organ-scale and drug-induced ion-channel blockage combination, we combine multiple machine learning techniques to create and train sex-specific arrhythmogenic drug risk classifiers. We first narrow the drug effect parameter space by studying the cellular drug sensitivity to various ion channel blockings. For this sensitivity analysis, we use multivariable logistic regression techniques, as they have been proven to be computationally more efficient than one-at-a-time parameter sensitivity analyses (Lee et al., 2013) and highly suitable for studying processes with binary outcomes (Morotti and Grandi, 2017). Next, we apply the developed sex-specific high-fidelity multiscale exposure-response simulators to quantify the risk of drug-induced arrhythmogenesis within the identified critical drug-induced ion channel blocking parameter space. To reduce the computational cost of exploring this parameter space, we develop and train multi-fidelity risk classifiers that have been shown to outperform single high-fidelity risk classifiers (Sahli Costabal et al., 2019b). More specifically, we combine low-fidelity mid-wall cardiomyocyte simulations and high-fidelity heart simulations to train a Gaussian process classifier that characterizes the probability of arrhythmogenicity based on the two most important ion channel blockage features for arrhythmogenesis. Using active learning, we maximize the information gained by each possible low-and high-fidelity sample we evaluate, keeping the computational costs of training our arrhythmogenic risk classifiers as low as possible.

Using our multi-fidelity arrhythmogenicity classification boundary, we estimate the arrhythmogenic risk of three drugs, a high, intermediate and low risk drug (Li et al., 2018), by computing the critical drug concentration at which arrhythmia will start developing. We select three drugs for which the concentration-block response curve is well-described (McMillan et al., 2017) for the two cardiac currents that have the most significant impact on arrhythmogenic risk prediction (section 2.3.1). The critical drug concentration is found at the intersection of the multi-fidelity arrhythmogenesis classification boundary and the two-dimensional concentration-block trajectory described by Equation (5). If the drug's concentration-block trajectory does not cross the risk boundary, that drug can be considered safe for the studied sex.

Figure 6 represents the male and female anti- or pro-arrhythmic sensitivity to drug-induced ion channel blocking. As can be seen in the upper two plots, the female midwall cells are more sensitive to drug-induced ion channel blocking than male midwall cells. For the same set of ion channel blocking samples B={βKr,βNa,βNaL,βCaL,βKs,βto,βK1}∈[0.0,0.95]10,000×7, we recorded 760 and 4,450 abnormalities for the male and female midwall cell, respectively.

The lower panel plots in Figure 6 depict the normalized marginal effects of drug-induced ion channel blocking on de- and repolarization abnormalities. For both male and female midwall cardiomyocytes, I_Kr blocking has a strong pro-arrhythmic effect, whilst I_CaL blocking has the largest anti-arrhythmic strength. The male and female normalized marginal effect of L-type Ca²⁺ channel blocking amounted to −0.272 and −0.190, respectively. As such, we selected drug-induced blocking of the rapid delayed rectifier potassium current, β_Kr, and the L-type Ca²⁺ current channel, β_CaL as the two main pro-arrhythmic and anti-arrhythmic input features β₊ and β₋ for our male and female arrhythmogenic risk classifier.

In training our multi-fidelity arrhythmogenic risk classifiers, we first trained a single fidelity de- and repolarization classifier for the male and female midwall cell, respectively. Figures A1, A2 summarize the de- and repolarization abnormality classification boundary delineation in the studied {β_CaL, β_Kr} parameter space. The upper panel plots showcase the initial exploration phase to train these Gaussian process classifiers. The lower panel plots depict the subsequent exploration and exploitation phase through active learning.

The subsequent training and development of the male multi-fidelity arrhythmogenic risk classifier is showcased in Figure 7. The upper panel plots showcase the evaluation of 10 high-fidelity evaluations of {β_CaL, β_Kr} on male drug-induced arrhythmogenesis. Left, the virtual electrocardiograms showcase the effect that various drug-induced ion channel blocking combinations have on the male heart. Here, only one exploratory sample (β_CaL = 6.7%, β_Kr = 92.6%) resulted in reentrant arrhythmia in the male heart. The other {β_CaL, β_Kr} combinations affected the QT interval, but did not lead to arrhythmogenesis. In the upper middle plot, the N_L = 50 low-fidelity evaluations are shown together with the first 10 exploratory high fidelity arrhythmogenicity classifications. On the upper right, the initial multi-fidelity Gaussian process risk classifier for male arrhythmogenesis is shown. Concomitantly, Figure A3 showcases a male single-fidelity risk classifier, only taking into account these high-fidelity arrhythmia development evaluations. Comparing the upper panel plots of Figure 7 with Figure A3, it can be seen that taking the low-fidelity classification data into account in training a Gaussian process classifier significantly aids the precision of the high-fidelity classifier with a limited amount of samples. This is the power of multi-fidelity Gaussian process classification. In the lower panel plots of Figure 7, we showcase the multiscale evaluation of 15 additional high-fidelity active learning samples. These samples allowed us to capture the bias in the low-fidelity predictions (see Equation 13) showcased in more detail in Figure A5 (left). The virtual electrocardiogram recordings of a subset of these active learning samples, four arrhythmic and four non-arrhythmic samples, are shown in the lower left and mid plots, respectively. Finally, the fully explored and exploited male multi-fidelity arrhythmogenic risk classifier is shown in the lower-right plot, with N_L = 50 and N_H = 25 low- and high-fidelity risk evaluations, respectively.

Figure 8 showcases the training of the female multi-fidelity arrhythmogenic risk classifier. In evaluating 10 high-fidelity Latin hypercube samples, five ion channel blocking samples drove the female heart to arrhythmogenesis, as shown in the upper middle plot. The electrocardiograms of two arrhythmic and two non-arrhythmic samples are shown in the upper left plots. The upper right plot depicts the initial female multi-fidelity drug-induced arrhythmogenicity classifier, taking into account all low-fidelity abnormality classification samples and the first 10 exploratory high-fidelity drug-induced arrhythmogenesis evaluations. Again, comparing this multi-fidelity classifier to the single-fidelity multiscale classifier shown in Figure A4 showcases the power of multi-fidelity Gaussian process risk classification. Through active learning, the low- to multi-fidelity bias (Equation 13) is inferred from 15 additional high-fidelity arrhythmogenic risk evaluations and depicted in Figure A5. Compared to the male heart, our results showcase a larger low- to multi-fidelity bias for the female heart. The lower left and mid plots in Figure 8 delineate the virtual electrocardiograms of four arrhythmogenic and four non-arrhythmogenic high-fidelity active learning sample evaluations, respectively. The final female multi-fidelity drug-induced arrhythmogenic risk classifier is shown in the lower-right plot.

Both male and female multi-fidelity drug-induced arrhythmogenicity classifiers highlight the pro-arrhythmic effect of I_Kr ion channel blocking and the anti-arrhythmic effect of I_CaL ion channel blocking. For the male heart, we predict drug-induced arrhythmogenicity at 81.7% I_Kr blocking when there is no I_CaL blocking. At 25.0, 50.0, and 75.0% I_CaL blocking, the critical β_Kr is 83.8, 88.6, and 94.6%, respectively. For 100% I_CaL blocking, no arrhythmia develops, regardless of the amount of I_Kr blocking. For the female heart, our risk classifier predicts drug-induced arrhythmogenesis at 43.7% blocking without any I_CaL blocking. For 25.0, 50.0, 75.0, and 100% I_CaL blocking, arrhythmia develops at 48.1, 52.2, 56.5, and 59.9% I_Kr blocking, respectively. Overall, the female heart can be expected to be significantly more prone to drug-induced arrhythmogenicity.

Figure 9 demonstrates how we use our male and female multifidelity arrhythmogenicity classification boundary to perform a sex-specific drug risk assessment. More specifically, we calculate the drug-induced arrhythmogenic risk for dofetilide, a high risk drug, chlorpromazine, an intermediate risk drug, and diltiazem, a low risk drug. For each of these drugs, the drug-specific color-coded block-concentration characteristics map onto a trajectory in the β_CaL/β_Kr plane. The intersection of this trajectory with our trained classification boundary defines the critical drug concentration at which arrhythmia can be expected to develop. For dofetilide, the block-concentration curve crosses the male and female arrhythmogenicity classification boundary at 26.0x and 3.5x the drug's effective free therapeutic plasma concentration, respectively. For chlorpromazine, the block-concentration curve does not cross the male classification boundary, signifying this drug can be considered safe for men. For women, the chlorpromazine block-concentration and classification boundary intersect at a critical concentration of 80.1x. For diltiazem, the block-concentration trajectory does not cross the male, nor the female arrhythmogenicity classification boundary, showcasing this drug to be safe for both sexes.

Until recently, the gold standard to assess pharmacological pro-arrhythmic risk consisted of assessing the potential of a drug (1) to cause pharmacological block of the rapid-delayed rectifier potassium I_Kr current encoded by the human ether-à-gogo related gene and (2) to prolong the QT interval in electrocardiographic animal and human studies. Although these biomarkers show good sensitivity, they are costly and have poor specificity, potentially blocking safe new drugs from ever reaching the market (Sager, 2008). In response to this problem, the Comprehensive in vitro Proarrhythmia initiative was launched (Sager et al., 2014). This incentive aimed to develop novel drug-induced arrhythmia biomarkers through a combined in vitro and in silico approach. in vitro, insights into the effect of drugs on multiple ion channels in the cardiomyocyte were collected. In silico, these insights were used to develop a mechanistic understanding how these ion channel blockings affect cardiac electrophysiology and function. Over the past decade, our in silico mechanistic understanding of the effects of drugs on cardiac electrophysiology has taken big leaps forward. As such, drug cardiotoxicity has been extensively studied in ventricular myocyte models (Mirams et al., 2011; Passini et al., 2017), transmural cable simulations (Moreno et al., 2013; Romero et al., 2018), planar and cubic tissue slabs (Kubo et al., 2017; Yang et al., 2020; Margara et al., 2021), and even ultra-high resolution, multiscale heart models (Wilhelms et al., 2012; Okada et al., 2015; Sahli Costabal et al., 2018a; Hwang et al., 2019). Within this paradigm, multiple groups have used such models to develop arrhythmogenic risk classifiers. These in silico augmented biomarkers showcase improved sensitivity and specificity with respect to the gold standard human ether-à-gogo and QT prolongation guidelines (Passini et al., 2017; Li et al., 2019). Currently, these in silico arrhythmogenicity biomarkers focus mainly on lower-fidelity isolated cardiac cell models (Lancaster and Sobie, 2016; Britton et al., 2017; Fogli Iseppe et al., 2021) or simplified cable simulations (Polak et al., 2018; Romero et al., 2018; Yang et al., 2020). The underlying motivation for such an approach is the role of cellular early afterdepolarizations and repolarization failures in providing a trigger for the development of arrhythmia. Nevertheless, arrhythmogenicity is not completely governed by, nor exclusively limited to, depolarization or repolarization abnormalities (Pugsley et al., 2015). Overall, the spatial dispersion of repolarization within the ventricular myocardium has been identified as the principal arrhythmogenic substrate (Antzelevitch and Burashnikov, 2011). A recent computational multiscale cardiac electrophysiology study showcased that the electrotonic coupling effect in tissue is an essential factor to predict drug effects on the living organ (Kubo et al., 2017). More specifically, computational 2D tissue slab results revealed no tachyarrhythmia in the presence of early afterdepolarizations at the cellular level, and arrhythmogenic induction in between the endocardial and midwall tissue layers, rather than in the midwall tissue itself. As such, it can be appreciated that accurate mechanistic understanding of arrhythmogenesis requires a high-fidelity multiscale modeling approach coupling the effect of drugs to subcellular ion channel activity, to cell-to-cell coupling at the tissue scale, and the tissue's three-dimensional heterogeneous and anisotropic organization at the organ scale.

The current lack of multiscale computional modeling in the development of novel in silico augmented arrhythmogenic risk classifiers can be associated to their computational cost. Whereas a single cell action potential takes seconds to compute on a single CPU, a coupled cell-tissue-organ scale exposure-response simulator can easily take multiple hours to compute on a high performance computing cluster (Towns et al., 2014; Sahli Costabal et al., 2018a). Upon developing a risk classifier that evaluates the arrhythmogenic susceptibility to a whole series of drugs at multiple drug concentrations, the computational burden of performing a multiscale evaluation for each case becomes too high. To overcome this limitation, our study took a different approach. Instead of evaluating the case-by-case drug- and concentration-specific response, we trained a risk classifier based on the most important drug-induced ion channel blockings. We used a combined multiscale modeling and machine learning approach, entailing logistic regression, single- and multi-fidelity Gaussian process classification, and active learning techniques. First, we quantified the principal role that drug-affected ion channel currents play in developing arrhythmia. Using these insights, we established a two-dimensional drug blocking parameter space in which we evaluated the arrhythmogenic susceptibility of various drug-induced ion channel blocking combinations. Next, we relied on Gaussian process classification to delineate the arrhythmogenicity border within the considered parameter space. Given the high computational cost of each multiscale evaluation, we set up a multi-fidelity framework. Here, we used cellular midwall cell evaluations as a low-fidelity proxy for arrhythmogenic risk classification. This low-fidelity classifier was subsequently used to inform the underlying shape of the multi-fidelity classifier. Doing so, we minimize the amount of high-fidelity evaluations within the studied parameter space, and still end up with a precise multiscale arrhythmogenic risk classification boundary. This multi-fidelity ion channel blocking classifier subsequently allows us to post-process the intrinsic arrhythmogenic risk for each possible drug and concentration at no additional computational cost, without losing accuracy of the underlying multiscale arrhythmic and non-arrhythmic classification. We took advantage of the probabilistic nature of our Gaussian process classifiers to implement an effective data acquisition via active learning strategies. These strategies sought a balance between parameter space exploration and classification boundary exploitation. Consequently, our methodology allowed us to maximized classifier accuracy under a constrained computational budget and provided a significant advantage over other classifiers including logistic regression and support vector machines.

About a century ago, sex differences in cardiac electrophysiology were reported for the first time (Bazett, 1920). Throughout the past two decades, these insights have matured into the recognition that female sex is an essential risk factor for multiple adverse cardiac events (Yarnoz and Curtis, 2008). Especially for drug-induced arrhythmogenesis, women turn out to be impacted twice as much as men (Makkar, 1993; James et al., 2007; Coker, 2008). Nevertheless, the effect of sex differences on cardiac electrophysiology and drug-induced arrhythmogenicity remain largely underexplored. With current sex-agnostic population-based models (Muszkiewicz et al., 2016; Li et al., 2019) being calibrated on in-vitro studies, which tend to be male-dominated (Ramirez et al., 2017), sex bias is expected to propagate through these novel in silico augmented arrhythmogenic risk classifiers. This study sought to take female sex into account as an independent biological variable by developing two sex-specific in silico augmented multiscale arrhythmogenic risk classifiers. To accomplish this, we first extended the multiscale envelope of studying sex-differences in cardiac electrophysiology beyond the cell or tissue level (Yang et al., 2017; Fogli Iseppe et al., 2021) up to the organ scale. Next, we used the developed framework to delineate male vs. female arrhythmogenic sensitivity to drugs.

Our male and female cell models were based on a high-throughput quantitative assessment of genome-scale sex differences in male and female human endo- and epicardial tissue (Gaborit et al., 2010). The resulting female endo- and epicardial action potential duration is significantly longer than the male action potential durations. Both the male and female endo- and epicardial action potential durations in this study are smaller than those computed in other studies (Yang and Clancy, 2012). Whereas other studies considered the baseline O'Hara-Rudy model and ion channel conductances to form the male baseline cell model (Yang and Clancy, 2012), our approach acknowledged the originally reported data population (O'Hara et al., 2011) and regarded the baseline model as a mixed 56% male / 44% female generalized model. Despite these differing views, our computed action potential durations fall well within the reported ranges based on experimental variability (Gaborit et al., 2010; Yang and Clancy, 2012). Similarly, the range of our sex-specific endo- and epicardial action potential durations are in agreement with reported populations of ventricular cell models (Britton et al., 2017). Our midwall cell action potential durations also fall within the same reported population variability. Averaged over the three cell types, our simulated female cells take 30% longer to repolarize than their male counterparts, which is consistent with the reported 29% relative female-to-male action potential differences for human ventricular myocytes (Verkerk et al., 2005). Focusing on the restitution behavior, our reported male and female action potential durations at 90% repolarization in Figure 3 agree favorably with previously reported experimental data for human tissue (Morgan et al., 1992; Drouin et al., 1995; Li et al., 1998; ten Tusscher et al., 2004; O'Hara et al., 2011).

At the organ scale, the combined effect of sex-specific cell-scale ion channel activity, tissue-scale conductivity and organ-scale geometry results in a shorter QRS and longer QT interval for women. Both results are in agreement with clinical population studies. Female vs. male QRS shortening of 5 ms, vs. 4 ms here, and a QT prolongation of 20 ms, vs. 63 ms here, have been reported in the literature (Vicente et al., 2014). The mismatch between a recorded 29% AP prolongation and a clinical QT prolongation of ‘only’ 2-6% has been hypothesized to be related to the mismatch between single isolated cell behavior and three-dimensional electrophysiologically coupled heterogeneous tissue (Verkerk et al., 2005). Indeed, our multiscale models showcase that a 30% action potential prolongation between both sexes at the cell scale only resulted in a male QT interval which was 15% shorter than the female QT interval at the organ scale. Nevertheless, this sex difference in QT interval duration is still on the higher end. This discrepancy seems to be related to our male multiscale heart model. The computed male QT interval of 348 ms corresponds to the 5th percentile of the clinically reported ranges for men (Asatryan et al., 2021), whereas our computed female QT interval of 411 ms aligns perfectly with the clinically reported range of 386–445 ms (Vicente et al., 2014). As our multiscale models demonstrated the dominant role that changes in ion channel activity have on the timing of the T wave end, there is a strong need for an in-depth experimental study on the sex-specific differences in functional ion channel activity of non-diseased human ventricular myocytes. Unfortunately, we are not aware of such data being currently available. Similarly, studies have shown that the inclusion of interventricular and apicobasal ion channel gradients at the tissue scale can further impact ECG morphology (Okada et al., 2011). Emerging electrocardiographic imaging techniques show great potential to study sex differences in healthy tissue-scale conductivity in more detail but remain challenging (Cluitmans et al., 2018; Andrews et al., 2019).

Given the high amount of ionic currents constituting the electrophysiological behavior of human ventricular cardiomyocytes (Equation 4), studying the drug-induced risk to develop arrhythmia requires the exploration of a large parameter space constituting different amounts of drug-induced blocking of each and every possible ion channel. To keep the parameter space computationally tractable, we focused on the seven most important ion channels for arrhythmogenic risk stratification, and used logistic regression to quantify their relative importance. The normalized marginal effects of drug-induced ion channel blocking on arrhythmic sensitivity in Figure 6 identify β_Kr and β_CaL as the key pro-arrhythmic and anti-arrhythmic ion channel blockings, respectively. This conclusion is consistent with previous sex-agnostic risk analyses (Crumb et al., 2016), and is thus found to hold true across men and women. Interestingly, our analysis highlights a relatively decreased protective role of L-type Ca²⁺ channel blocking in women. The higher amount of recorded de- and repolarization abnormalities confirmed the higher susceptibility of female cardiomyocytes to drug-induced arrhythmogenicity. These results agree well with experimental exploratory studies on cell-scale sex differences in drug-induced arrhythmogenicity (Liu et al., 1999; Verkerk et al., 2005). The male and female multi-fidelity arrhythmogenic risk classifiers in Figures 7, 8, respectively quantify this differing risk with increased fidelity, as shown in Figure A5. Overall, we found the female heart to demonstrate arrhythmogenicity at lower drug-induced I_Kr and I_CaL ion channel blocking than the male heart. Interestingly, our previous work on sex-agnostic arrhythmia risk assessment in the heart showcased an arrhythmogenic risk classification boundary in between the male and female arrhythmia risk classification boundary developed in this study (Sahli Costabal et al., 2019c). As such, we conclude that a generalized sex-agnostic arrhythmia risk classification underestimates and overestimates the cardiac toxicity of drugs for women and men, respectively. This directly puts women at higher risk for drug-induced arrhythmogenicity events, explaining the higher incidence reported in women (Makkar, 1993; James et al., 2007; Coker, 2008).

In applying our novel sex-specific arrhythmogenic risk classifier to a high, intermediate and low risk drug, we quantify this increased risk for women in more detail. For dofetilide, a class III anti-arrhythmic agent, both the male and female arrhythmia risk classifier confirm the general notion that dofetilide can have dramatic consequences if not dosed correctly (Briceño and Supple, 2017). For men, our risk classifier predicts a spontaneous transition from a sharp but smoothly propagating excitation pattern into rapid, irregular, asynchronous activation patterns at a critical concentration of 26.0x. For women, the same risk is predicted at 3.5 times the drug's free therapeutic plasma concentration. These results agree well with clinical trials where female sex was associated with three-fold higher odds of dofetilide discontinuations or dose reductions relative to the male sex (Pokorney et al., 2018). Most dosage reductions led to half of the recommended dosage for women. Interestingly, women were highly underrepresented in original clinical trials assessing the safety of dofetilide, only accounting for 28 and 16% of the total amount of enrolled patients (Køber et al., 2000; Singh et al., 2000). For chlorpromazine, an antipsychotic drug, our female arrhythmogenic risk classifier estimated a risk for arrhythmogenesis at 80.1x concentration, whilst for men no arrhythmogenicity was predicted. As expected from such a high critical risk concentration, chlorpromazine-induced arrhythmogenicity can be expected to be uncommon. Indeed, a comprehensive literature search spanning four decades of clinical case report data identified only seven published cases of chlorpromazine-associated ventricular arrhythmia. All these cases involved women (Hoehns et al., 2001). Finally, for diltiazem, a calcium channel blocker used to manage blood pressure and chest pain, the drug's concentration-block trajectory does not cross our male nor female multi-fidelity arrhythmogenic risk classification boundary. Consequently, we predict no arrhythmogenesis for diltiazem and consider this drug to be safe, both for men and women. This risk assessment corresponds well with diltiazem's ‘low/no arrhythmia risk’ classification by a team of clinical cardiologists and electrophysiologists based on publicly available data and expert opinion (Colatsky et al., 2016). Additionally, the arrhythmogenic safety of diltiazem was also confirmed by recent sex-agnostic population-based arrhythmia risk classifiers from other research groups (Lancaster and Sobie, 2016; Li et al., 2019). Importantly, these sex-specific drug-induced arrhythmogenic risk assessments assume no other medications to be taken concomitantly with these drugs.

In this study, we first build a multiscale mechanistic understanding of arrhythmogenesis in the male and female heart, and subsequently use these computational insights to evaluate the sex-specific drug-induced arrhythmogenic risk. This approach is inherently different from a recent first approach toward sex-specific drug-induced arrhythmogenicity classification (Fogli Iseppe et al., 2021). In this study, the authors focus on the in silico computed effect of drugs on male and female human epicardial cardiomyocytes. Using statistical learning techniques, they identify the key synthetic action potential biomarkers contributing to the most accurate prediction of arrhythmogenicity outcomes for men and women specifically. This approach relied on a ground truth classification assumption that risky drugs are dangerous for men and women, and safe drugs are safe for both men and women. These classifications were deduced from reported adverse event analyses performed within the Adverse Drug Event Causality Analysis (Woosley et al., 2017). With female sex reported to be historically highly underrepresented in clinical studies (Vitale et al., 2017), such an assumption could potentially be problematic, especially for older drugs. Our study offers the benefit of using a full multiscale framework to inform arrhythmogenic risk from a mechanistic understanding. Apart from this differing approach to arrhythmogenic risk classifier development, our study also takes into account the effect of midwall cells in an individual's predisposition to arrhythmogenesis (Drouin et al., 1995; O'Hara et al., 2011) and did not assume the baseline O'Hara-Rudy model to represent a purely male endocardial cell model, as discussed in section 2.2.1. Consequently, these differing approaches render a one-on-one comparison between our studies difficult. Nonetheless, for chlorpromazine, the only drug that was studied in both studies, both our studies classify this drug safe for men and at medium risk for women. Based on our sex-specific arrhythmia risk classifiers in which the male arrhythmogenic {β_CaL, β_Kr} risk zone ⊂ the female arrhythmogenic {β_CaL, β_Kr} risk zone, our classifiers do not predict any drug to have a higher risk for women than for men. Interestingly though, the alternative approach identified specific drugs that are safer for women than for men (Fogli Iseppe et al., 2021). This disagreement might be associated with the current ambiguity on functional sex differences in Ca²⁺ handling (Parks and Howlett, 2013; Parks et al., 2014) which led the authors to disregard the genomic sex differences in sodium potassium I_NaK pump and Ca²⁺ uptake channel activities and relatively upscale the female sodium calcium exchange currents we deduced in Table 1. Currently, experimental data on Ca²⁺ in healthy human myocardium are lacking, and further investigation on these functional sex differences is warranted to improve our sex-specific arrhythmogenic risk classifiers in the future. Additionally, limiting our risk classifier to only take β_CaL and β_Kr into account might not uncover this behavior. It can for example be seen from Figure 6 that β_NaL has a stronger anti-arrhythmic effect for female midwall cardiomyocytes than it does for men, and thus including this feature as a third drug-induced ion channel blocking input feature to our arrhythmia risk classifier might lead to male and female three-dimensional arrhythmogenic risk zones that do not completely overlap. Apart from these discrepancies, the overall conclusion is the same: including sex as a new independent factor in preclinical cardiotoxicity risk assessment is crucial to avoid potentially life-threatening consequences for the female population (Chorin et al., 2017). Given the absence of reliable large-scale arrhythmogenic risk assessments for women specifically, and the male dominance in clinical studies, our study forms an important first step toward mechanistically uncovering the role that sex differences on the subcellular, cellular, tissue, and organ scale play in drug-induced arrhythmogenicity. An improved understanding of these sex-specific mechanisms will be crucial to provide new therapeutic approaches that do no longer put women at increased risk.

Although our proposed methodology holds great promise to rapidly assess the sex-specific risk of a new drug, without relying on clinically reported adverse event occurrence, it has a few limitations: first, our sex-specific multiscale exposure-response simulators are only as good as their input. In the long term, more sex-specific human cell and tissue experiments are needed to fine-tune the cell- and tissue-scale sex differences in ion channel activity and conductivity currently deduced from genomic data. Such experimental data would also be highly desirable to resolve the current debate on the existence or non-existence of sex differences in Ca²⁺ handling. The lacking human experimental data for more in-depth sex-specific validation of our multiscale simulation outcomes suggest important avenues for further studies. Novel developments of male and female hIPSC-derived cell lines might provide an interesting route to study this further (Huo et al., 2019). Second, our developed risk classifiers currently focused on the risk of drug-induced I_Kr and I_CaL blocking. Even though we identified these ion channel blockings to be the most critical channels for drug-induced arrhythmogenesis for both the male and female heart, arrythmogenic risk stratification for drugs that mainly target other channels might require including additional ion channel blockings to our risk classifiers. As our results in Figure 6 showcase, extending both the male and female risk classifiers to take into account drug-induced I_Ks and I_NaL blocking would be the most logical next step. Third, given the role that the excitation rate has on the electrophysiological behavior of the human heart, we aim to extend our classifiers to take into account heart rate in our future work. Fourth, at this point, we did not take into account male and female population variability. As additional experimental data becomes available, more in-depth sex-specific population studies form an interesting next step. We have developed efficient frameworks to quantify and propagate such uncertainty through computational models in the past (Peirlinck et al., 2019; Sahli Costabal et al., 2019a), and look forward to apply these techniques to this problem. Prior to this, a critical and logical next step would be to validate our method using our own independent experiments with human adult cardiomyocytes, and ideally, healthy human volunteers. Ultimately, with a view toward precision cardiology, this sex-specific approach forms an important initial step toward identifying the optimal course of care for each individual patient based on personalized block-concentration characteristics and personalized cardiac heart models (Trayanova, 2018; Peirlinck et al., 2021; Rodero et al., 2021).