Healthy human cardiomyocytes (hCMs) were modeled using the O’Hara et al. endocardial ventricular model (O'Hara et al., 2011) with sodium current dynamics modified as detailed by Mora et al. (Mora et al., 2018). To model non-fibrotic HF cardiomyocytes, several parameters were modified as detailed elsewhere by Mora et al. (Mora et al., 2018) to reflect the underlying pathophysiology of HF with reduced ejection fraction. HF in fibrotic myocardium was modeled as coupling of an individual HF with reduced ejection fraction cardiomyocyte to 5 fibroblasts with a gap junctional conductance of 1 nS, as detailed in our previous work (Phillips et al., 2021) (adapted from Mora et al. (Mora et al., 2018)). Pacing frequency was set at 1 Hz, and all simulations were run for 500 beats to achieve a steady state solution, in accordance with our previous work (Phillips et al., 2021).

Theoretical candidate non-excitable custom cells were developed for the correction of failing cardiomyocytes via heterocellular coupling, expressing membrane channels from hCICs, hMSCs, and cardiac fibroblasts (CFs).

Ionic currents from hCIC channels—including a sodium current, inward rectifying potassium current, large-conductance calcium-activated potassium current, and transient outward current—were modeled as described in our recent work (Phillips et al., 2021). Ionic currents from hMSC channels—including a transient outward current, sodium current, L-type calcium current, calcium-activated potassium current, and delayed rectifier potassium current—were modeled as described in our previous work (Mayourian et al., 2016). Finally, ion currents from CF channels—including the delayed rectifier potassium current, inward rectifier potassium current, and sodium/potassium pump current—were modeled as described elsewhere (MacCannell et al., 2007). Total ionic flow through the custom cell was subsequently calculated as the sum of the ionic flows through hCIC, hMSC, and CF channels, where each channel’s maximal conductance was scaled as described below as per our genetic algorithm.

The capacitance of all custom cells was set to 40pF, to reflect the approximate capacitance of hMSCs (Mayourian et al., 2016) and hCICs (Phillips et al., 2021), and to reflect their smaller size compared to hCMs (O'Hara et al., 2011). During simulations, intracellular ion concentrations of all custom cells were assumed to be held constant.

Another type of theoretical candidate custom cell was developed that expressed healthy hCM membrane channels in addition to the previously described non-excitable ion channels.

Ionic currents from healthy hCM channels—including the fast sodium current, late sodium current, transient outward current, L-type calcium current, rapid delayed rectifier potassium current, slow delayed rectifier potassium current, and inward rectifier potassium current—were modeled as described in Mora’s modification (Mora et al., 2018) of the O’Hara et al. human endocardial ventricular cardiomyocyte model (O'Hara et al., 2011). Total ionic flow through the custom cell was subsequently calculated as the sum of ionic flows through hCIC, hMSC, CF, and hCM channels, where each channel’s maximal conductance was scaled as described below as per our genetic algorithm.

Similar to above, the capacitance of all custom cells was set to 40pF, and the intracellular ion concentrations were held constant.

A human cardiomyocyte (non-fibrotic HF, fibrotic HF, or healthy) was coupled to custom cells using the following pair of equations, similar to our previous work (Phillips et al., 2021) and as described elsewhere (Xie et al., 2009): − C C M d V C M d t = I i o n ,   C M + n c u s t o m G g a p V C M − V c u s t o m (1) − C c u s t o m d V c u s t o m d t = I i o n ,   c u s t o m + G g a p V c u s t o m − V C M (2) where: C_CM and C_custom represent the membrane capacitance of the cardiomyocyte and custom cell, respectively; V_CM and V_custom represent the membrane voltage of the cardiomyocyte and custom cell, respectively; t is time; I_ion,CM and I_ion,custom represent total ionic flow through the cardiomyocyte and custom cell, respectively; n_custom is the integer number of coupled custom cells per cardiomyocyte; and G_gap is gap junctional conductance between the cells. Custom cells were only coupled to cardiomyocytes (i.e., custom cells were not also coupled to fibroblasts).

A genetic algorithm for parallel computational model fitting was implemented as described elsewhere (Groenendaal et al., 2015; Devenyi et al., 2017). An initial population of 2500 custom cell models was generated for each in silico experiment with varying ion channel maximal conductances. Maximal conductance parameters for each ion channel were defined as the original value multiplied by a scaling factor obtained by pseudo-random sampling from a uniform distribution ranging from 0.01% to 1000%. Gap junction conductance between the custom cell and the failing cardiomyocyte was similarly obtained by pseudo-random sampling from a uniform distribution ranging from 0.01% to 1000% of its designated baseline value of 1 nS. The upper limit for gap junction conductance was selected to minimize computational cost, as our recent work showed that when coupling cardiomyocytes to non-excitable cells, gap junctional conductance values greater than 10 nS had minimal additional effect on AP and CaT waveforms (Phillips et al., 2021). The scaled number of heterocellular coupled custom cells per cardiomyocyte was pseudo-randomly sampled from a discrete uniform distribution ranging from 1 to 5 (baseline value of one heterocellular coupled custom cell per cardiomyocyte). The upper limit of 5 was selected as a practical consideration of therapeutic delivery. For each custom cell model, the AP and CaT waveforms were extracted from the final (i.e., 500th) simulated beat.

Subsequent to model creation, the genetic algorithm was run for 5 rounds, or “generations”, aiming to achieve at least 10% acceptance of custom cells. Further rounds were purposely not pursued, as the solution space would further converge (whereas a wide solution space was of interest in this study). Each generation of the genetic algorithm consisted of 4 steps: selection, crossover, mutation, and elitism (Bot et al., 2012). During selection, models were randomly paired, and the AP and CaT root mean square errors (RMSE_AP and RMSE_CaT, respectively) between the custom cell-treated HF cardiomyocyte and a healthy cardiomyocyte for the first 500 ms of the final beat were defined as the objective functions. All RMSEs were normalized to the RMSE between a healthy and untreated failing cardiomyocyte. The model achieving a greater number of objective function values less than 50% (where RMSE_AP is the first objective function and RMSE_CaT is the second objective function) was kept. In the case of a tie, the model with the lower overall percent RMSE between the two waveforms was chosen. Selection was completed twice per generation such that the resulting number of models advancing remained at 2500, with each original model advancing either 0, 1, or 2 times post-selection. After selection, models that advanced were randomly paired once more for crossover. Each pair had a 90% chance of undergoing crossover, in which each parameter value then had a 50% chance of crossing over or exchanging between models. After crossover, each model parameter had a 1% chance of undergoing random mutation, where it would be scaled by a log-normal distribution with a standard deviation of 20% such that >99.9% of mutations would fall between a 0.5-fold and 2-fold change. Finally, an elitism strategy was employed such that the top 10% of models from each generation replaced the bottom 10%.

Following genetic algorithm adjustment, accepted models from the final generation were determined based on AP and CaT as compared to healthy cardiomyocytes. Specifically, models in the final generation were accepted if RMSE_AP and RMSE_CaT in the first 500 ms of the final beat were both less than 50% of their respective RMSEs between a healthy and untreated HF cardiomyocyte.

Following genetic algorithm adjustment, several metrics for electrophysiology and calcium handling in accepted models were calculated, including action potential duration at 50% (APD₅₀) and 90% (APD₉₀) repolarization, peak membrane voltage, maximum upstroke velocity, resting membrane potential, maximum systolic calcium concentration, maximum diastolic calcium concentration, calcium transient amplitude, and calcium relaxation time constant at 50% (τ₅₀) and 90% (τ₉₀) decay. These metrics were used to further interrogate the effects of custom cell coupling. Note statistical analyses were deferred given the nature of this work whereby arbitrarily large sample sizes could be achieved, which is in accordance with our previous studies (Mayourian et al., 2016).

First, to investigate whether custom cells with non-excitable cell channels can restore a healthy hCM phenotype, we created a population of 2500 custom cell models that variably express ion channels from hCICs, hMSCs, and CFs. The genetic algorithm was subsequently used to promote selection of models that restore failing electrical and calcium cycling behavior back towards a healthy phenotype.

In non-fibrotic HF cardiomyocytes, heterocellular coupling with the population of custom cells led to no acceptable solutions (Figure 1A). Similarly, restoration was unable to be achieved in fibrotic HF cardiomyocytes (Figure 1B). While custom cell coupling was able to correct AP waveforms in a majority (90%, or 2245/2500) of non-fibrotic HF models in the genetic algorithm final generation, it was unable to correct AP waveforms (0%, or 0/2500) in the fibrotic HF models. CaT was not restored in either HF condition.

Next, we tested whether custom cells expressing channels from cardiomyocytes can restore a healthy phenotype. To do so, we created a population of 2500 custom cell models that variably express ion channels from healthy hCMs, in addition to hCICs, hMSCs, and CFs (Figure 2). The genetic algorithm was subsequently used to promote selection of models that restore electrical and calcium cycling behavior back towards a healthy phenotype.

The addition of excitable cell channels to custom cells was insufficient to restore both a healthy AP and CaT cardiac phenotype in non-fibrotic HF cardiomyocytes (Figure 2A). While 1963/2500 (79%) custom cells met acceptance criteria for the AP, only 4/2500 (0.2%) of the custom cells had <50% RMSE_CaT (Figure 2A), whereby the CaT peak was restored by as much as 55% from non-fibrotic HF to healthy conditions.

In contrast, coupling a population of custom cells to fibrotic HF cardiomyocytes led to solutions that satisfied acceptance criteria for both AP and CaT within Generation 0 (Figure 2B). By Generation 5, 266/2500 (11%) custom cells met the acceptance criteria (Figure 2B). Notably, the AP morphology of accepted solutions differed from healthy waveforms, whereby the upstroke velocity and peak voltage appeared lower, phase 2 was flatter, and phase 3 repolarization was more rapid (Figure 2B). In addition, the accepted CaT appeared to have lower calcium amplitude, higher diastolic calcium, and slower relaxation than healthy hCMs (Figure 2B). Figure 3 corroborates these findings and quantifies the differences in standard AP and CaT waveform metrics for treated HF fibrotic cardiomyocytes in comparison to healthy and untreated fibrotic HF hCMs.

Distributions of accepted custom cell parameters are shown in Figure 2C. Accepted solutions were dependent on gap junction conductance and number of coupled custom cells parameters, which spanned a narrow range with median values of 7.2 [IQR 6.3-8.3] and 5 [IQR 4-5], respectively. While a majority of ion channel maximal conductance parameters spanned the entire search space, hCM L-type calcium current (I_CaL) and inward rectifier K⁺ current (I_K1) spanned a narrow range in order to correct the fibrotic HF phenotype with a median maximal conductance scaling of 4.3 [IQR 3.6-5.0] and 3.0 [IQR 2.1-3.8], respectively (Figure 2C).

Next, we sought to understand the key properties of accepted custom cells that led to restoration of failing fibrotic HF cardiomyocytes in Figure 2. Based on the narrow solution space for select parameters in Figure 2C, we hypothesized I_CaL and I_K1 were critical for phenotype restoration of failing cardiomyocytes. To test this, we repeated the genetic algorithm from Figure 2, but without I_CaL and I_K1 as candidate channels for the custom cell (Supplementary Figure S1). Indeed, in this case, no custom cells met the acceptance criteria (Supplementary Figure S1), suggesting either or both of these select channels (i.e., I_CaL and/or I_K1) are critical to creating custom cells that can effectively correct the fibrotic HF phenotype. We next created a population of 2500 custom cells including only the following parameters of interest: I_CaL, I_K1, gap junction conductance, and number of coupled custom cells.

As shown in Figure 4, coupling custom cells that incorporate only these 4 parameters to fibrotic HF cardiomyocytes satisfied acceptance criteria for both AP and CaT in Generation 0 (Figure 4A). By Generation 5, 579/2500 (23%) custom cells met acceptance criteria (Figure 4A). The AP morphology of accepted solutions more closely resembled healthy waveforms in comparison to custom cells from Figure 2B where all channels were candidates; this is corroborated by the RMSE, where Figure 2B custom cells with all channels as candidates had a median RMSE_AP of 0.45 [IQR 0.43-0.47] and median RMSE_CaT of 0.47 [IQR 0.45-0.49], in comparison to Figure 4A custom cells with select channels (i.e., I_K1 and I_CaL) which had a median RMSE_AP of 0.37 [IQR 0.32-0.41] and median RMSE_CaT of 0.43 [IQR 0.40-0.45] (Supplementary Figure S2). While the acceptance bounds for I_K1, gap junction conductance, and number of coupled custom cells were nearly unchanged when considering select channels (Figure 4B) rather than all channels (Figure 2C), interestingly, I_CaL had a lower and more restricted acceptance range (Figure 4B); this could be explained by I_CaL not having to counteract a wide range of repolarization currents which were previously included as candidate channels.

We next explored if incorporating only I_CaL, or only I_K1, would be sufficient to create accepted custom cells. As shown in Figure 4C, when incorporating only I_K1 as a candidate channel, accepted custom cells closely resembled incorporating both I_CaL and I_K1 with a similar RMSE (Supplementary Figure S2) and acceptance parameter range (Figure 4D). In comparison, when incorporating only I_CaL as a candidate channel (Figure 4E), accepted custom cells had a larger RMSE (Supplementary Figure S2) despite a similar acceptance bound (Figure 4F). The corresponding raw accepted parameter solutions for ion channels and their relevant experimental metrics are shown in Supplementary Table S1. Interestingly, the I_CaL-only custom cell solution had minimal I_CaL activity, and appeared to resemble a passive-cell type; however, running the genetic algorithm on passive cells alone led to no accepted solutions (Supplementary Figure S3). Altogether, these findings suggest cardiomyocyte inward rectifier K⁺ channel is most critical to accomplish the phenotypic rescue.

We next calculated single-cell AP and CaT metrics associated with arrhythmogenesis to explore potential adverse risks in our accepted custom cells incorporating I_CaL-only, I_K1-only, or both. The similarities between incorporating I_K1 only versus both I_K1 and I_CaL into custom cells were further confirmed by these metrics (Figure 5). In contrast, as shown in Figure 5, the custom cell only incorporating I_CaL had a longer APD₉₀ and increased APD triangulation (i.e., APD₉₀ minus APD₅₀) compared to I_K1-only and I_K1+I_CaL custom cells, which are surrogate markers of pro-arrhythmia. All custom cells had similarly decreased upstroke velocity, a surrogate marker of conduction velocity. Finally, I_CaL had a depolarized resting membrane potential compared to I_K1-only and I_K1+I_CaL custom cells, also a surrogate marker of pro-arrhythmia. Altogether, these single-cell metrics suggest custom cells incorporating only I_CaL may pose a greater risk of an arrhythmogenic substrate compared to I_K1-only and I_K1+I_CaL custom cells, further supporting the incorporation of cardiomyocyte I_K1 only in custom cells.

Next, we aimed to gain insight into the mechanisms by which our accepted custom cells were theoretically able to restore healthy AP and CaT waveforms in fibrotic HF cardiomyocytes. To do so, we first explored individual cardiomyocyte ionic currents during the action potential when coupled to 25 randomly selected accepted custom cell solutions (Figure 6).

As shown in Figure 6A, inward rectifier activity of failing cardiomyocytes when coupled to accepted custom cells were more representative of healthy cardiomyocytes; delayed rectifier activity was intermediate to untreated fibrotic HF cardiomyocytes and healthy cardiomyocytes (Figures 6B, C). Interestingly, as shown in Figure 6D, L-type calcium current was more robust in fibrotic HF cells treated with accepted custom cells. Given that this occurred as early as phase 0 of the action potential, this could likely be explained by achieving a favorable resting membrane potential to maximize channel kinetics. In addition, the later activity of the L-type calcium channel is intermediate to fibrotic HF cardiomyocytes and healthy cardiomyocytes, likely in the setting of prolonging phase 2 (Figure 6D inset). By increasing L-type calcium channel activity altogether, this is presumably increasing calcium-induced calcium release (with increased ryanodine receptor activity in phase 0 as shown in Figure 6F). This is further corroborated by increased sodium-calcium exchanger activity during repolarization indicative of more intracellular calcium being extracted (Figure 6E).

To gain further insight, we also explored the membrane voltage and gap junctional current for the same randomly selected accepted custom cells during the cardiomyocyte action potential (Supplementary Figure S4). As shown in Supplementary Figure S4C, the accepted custom cell acts predominantly as an electrical sink during phases 0, 1, late 3, and 4. In doing so, the peak and plateau membrane voltages of the cardiomyocyte are reduced (Supplementary Figure S4A), which may be more favorable for achieving larger magnitude L-type calcium channel currents based on the established current-voltage relationship (O'Hara et al., 2011), which may in turn also contribute to improved calcium handling.

Finally, given that custom cells expressing I_K1 alone showed promise for correcting fibrotic HF cardiomyocyte phenotype with minimized RMSE and AP single-cell markers of arrhythmogenesis, we explored the translational potential of genetically engineering standard cell-based cardiotherapeutics—namely hMSCs and hCICs—to express this cardiomyocyte channel.

To test this, we started with either hMSCs or hCICs, and varied the following parameter variables: number of cells, gap junction conductance, and I_K1 maximal conductance. As shown in Figure 7, while no candidate cells were accepted, the minimum RMSE_AP of hMSCs decreased from 81% to 57% when expressing I_K1, and the minimum RMSE_CaT decreased from 79% to 43% when expressing I_K1 (Figures 7A, B; pink waveforms). Similarly, in hCICs, the minimum RMSE_AP decreased from 87% to 54% when expressing I_K1, and the minimum RMSE_CaT decreased from 79% to 38% when expressing I_K1 (Figures 7C, D; pink waveforms). The corresponding raw ion channel parameter values to achieve minimum RMSE and their relevant experimental metrics are shown in Supplementary Table S1.

We further explored the effects of these custom engineered hMSCs and hCICs achieving minimum total RMSE (see Figure 7 pink waveforms) on single-cell metrics associated with arrhythmogenesis—namely APD, resting membrane potential, and upstroke velocity. As shown in Figure 8, hMSCs and hCICs expressing I_K1 had APD₉₀, APD₅₀, and therefore APD triangulation closer to healthy conditions (Figure 8A). In addition, resting membrane potential better resembled healthy conditions (Figure 8C). Upstroke velocity for hCICs containing I_K1 was increased compared to hCICs without I_K1; upstroke velocity of hMSCs with versus without I_K1 was nearly unchanged (Figure 8B).

The concept of delivering somatic cells genetically engineered to exert a specific electrical influence on host cardiomyocytes has been previously hypothesized and experimentally investigated (Yankelson et al., 2008; Kirkton and Bursac, 2012; McSpadden et al., 2012; Nussinovitch et al., 2014; Kostecki et al., 2021). Several proof-of-concept studies have demonstrated the capability of engineered cells to alter cardiomyocyte AP, CaT, and arrhythmogenesis. For example, Yankelson et al. showed that mouse fibroblasts genetically engineered to express inward or delayed rectifier potassium channels hyperpolarized resting membrane potential, reduced action potential amplitude, and prolonged the cardiac effective refractory period when delivered to pig myocardium (Yankelson et al., 2008). The Bursac group later showed that coupling genetically engineered human embryonic kidney cells with variable inward rectifier and gap junction channel expression led to distinct effects on cardiac AP shape and initiation (McSpadden et al., 2012). These findings motivated more recent work by Kirton and Bursac, where delivering genetically engineered human embryonic kidney cells expressing inward rectifier, sodium, and gap junction channels in an in vitro model of fibrotic neonatal rat myocardium led to improved twitch force amplitude (Kirkton and Bursac, 2012).

These experimental findings complement our computer simulations. For example, interestingly our genetic optimization algorithm also identified cardiomyocyte inward rectifier channel as critical to improving AP and CaT waveforms in failing hCMs. The expression of these channels in custom cells also led to hyperpolarized resting membrane potential, reduced action potential amplitude, and higher calcium transient amplitude (a surrogate of twitch force).

Our findings provide theoretical insight into approaches to improve translational cell-based cardiac therapies. First, we demonstrate that custom cells that express cardiomyocyte ion channels are capable of correcting failing cardiomyocyte behavior, and that expression of cardiomyocyte inward rectifier and gap junction connexin-43 activity may be sufficient to accomplish this. These findings complement ongoing investigations; for example, the Bursac group has recently created a pipeline to generate and customize cell-based therapies (Nguyen et al., 2018), and has demonstrated that primary human fibroblasts expressing sodium channels, an inward-rectifying potassium channel, and connexin-43 can rescue conduction slowing in an in vitro model of cardiac interstitial fibrosis (Nguyen et al., 2016). However, in contrast to these studies, our custom cells decreased upstroke velocity (as our acceptance criteria were RMSE of AP and CaT, leading to phase 0/upstroke velocity not being prioritized in the genetic algorithm). This may be mitigated if custom cells also incorporated sodium channels (like the aforementioned experimental studies), a potential area of future investigation.

Second, our findings whereby custom engineered hMSCs and hCICs expressing inward rectifier channels reduced the RMSE in failing cardiomyocytes provide an avenue for improving current cell-based therapeutics. Predominant hMSC and hCIC cardiotherapeutic mechanisms of action involve paracrine signaling and heterocellular coupling (Menasche, 2018). As shown in our recent simulations, heterocellular coupling alone in standard hMSCs and hCICs was insufficient to correct a failing hCM phenotype (Phillips et al., 2021), motivating the current effort to design enhanced cell-based therapeutics. Our simulations suggest that incorporation of an inward rectifier channel improves a failing cardiac phenotype, and may minimize arrhythmogenic risk based on our single-cell metrics; these findings motivate future efforts to experimentally test custom hMSCs, hCICs, or CardioChimeras expressing an inward rectifier channel to optimize heterocellular coupling effects of delivered cells. Notably, such heterocellular coupling should be considered in the context of established paracrine effects of these therapeutic cells (Mayourian et al., 2017; Mayourian et al., 2018b).

Third, our computational study provides a framework for optimizing other therapeutic cell types. For example, in vivo delivery of human induced pluripotent stem cell-derived cardiomyocytes to the failing heart remains an active area of investigation (Dhahri et al., 2022). Limitations of this modality include an immature cardiac phenotype and poor electrical integration with host tissue (Dhahri et al., 2022). Our computational framework can be readily modified to optimize ion channel expression to improve the therapeutic potential of such stem cell-derived hCMs.

In addition to providing avenues of future translational work, our simulations provide insight into healthy and failing cardiomyocyte electrophysiology. First, our simulations suggest L-type calcium current was more robust in fibrotic HF cells treated with accepted custom cells; this could be explained by achieving a favorable resting membrane potential to maximize channel kinetics, or by achieving a more favorable peak cardiomyocyte membrane voltage that would correspond to larger magnitude currents based on the L-type calcium channel current-voltage relationship. By increasing L-type calcium channel activity, this would presumably increase calcium-induced calcium release.

Interestingly, our candidate custom cells expressing all non-excitable and cardiomyocyte ion channels led to a larger RMSE than expressing only cardiomyocyte inward rectifier and L-type calcium channel. In the latter experiment, AP phase 1-3 outward currents are set to zero in custom cells, leading to no competing activity against inward L-type calcium channel activity; by doing so, we reveal the minimal role of L-type calcium channel needed in our accepted custom cells. Interestingly, to our surprise, nearly passive cells with minimal L-type calcium channel activity are also accepted, and passive cells alone are nearly accepted, albeit with adverse prolonged APD triangulation and APD₉₀ at risk for arrhythmogenesis.

In this study, the wide range of candidate ion channels presented challenges for the unique and robust identification of optimal channels for expression in custom cells. Herein, we demonstrate the effectiveness of a genetic algorithm-based method to design custom cells that maximize HF correction. Inspired by Darwinian evolution and natural selection, genetic algorithms are an evolutionary metaheuristic method executed in a stochastic combinatorial manner (Akwaboah et al., 2021). One advantage of this method (compared to gradient-based optimization methods) is that it is less susceptible to getting stuck at local minima, often converging at the globally optimal solution even with high-dimensional parameter spaces (Gulsen et al., 1995). Similar to our approach, genetic algorithms have previously been applied to cardiac electrophysiology models using simple perturbations of maximum channel conductance to generate a desired AP morphology (Syed et al., 2005; Bot et al., 2012; Tomek et al., 2019; Smirnov et al., 2020).

To ensure biologically relevant solutions, we imposed physiologic bounds on the searchable parameter space, guided by following previous work (specifically for gap junction conductance) (Valiunas et al., 2004; Phillips et al., 2021) or by scaling established ion channel maximal conductances for non-excitable cells (Phillips et al., 2021) or cardiomyocytes (O'Hara et al., 2011; Mora et al., 2018). In addition, accepted solutions were restricted to no more than 5 coupled cells per cardiomyocyte, which would be challenging to exceed in clinical practice (although theoretically it is within the modeling capabilities).

This study has some recognized limitations. First, our candidate channels were limited to previously modeled non-excitable and cardiomyocyte channels. To our knowledge, we used the only published models for hMSCs (Mayourian et al., 2016) and hCICs (Phillips et al., 2021). Other cardiomyocyte and fibroblast models are available; however, to our knowledge only the combination of the O’Hara (O'Hara et al., 2011) and MacCannell (MacCannell et al., 2007) models have been used to develop a human fibrotic HF cardiomyocyte model, prompting us to use these models. We note that hMSC, hCIC, and CF leakage currents were excluded in our genetic algorithm, as they are often mathematical constructs to match electrophysiological measurements rather than biologically relevant channels, and by study design we aimed to express physiologic channels.

Second, by experimental design, in this study we solely focused on heterocellular coupling effects of custom cells, including engineered hCICs and hMSCs. However, it is well established that hMSC and hCIC cardiotherapeutic mechanisms of action also involve paracrine signaling (Menasche, 2018). For example, previous work showed that hMSC paracrine signaling alters excitation-contraction coupling of cardiomyocytes, causing increased expression of key calcium cycling genes—including sarcoendoplasmic reticulum Ca²⁺-ATPase and L-type calcium channels—and thereby increasing cardiomyocyte calcium transient amplitude and contractility (DeSantiago et al., 2012; Mayourian et al., 2017; Mayourian et al., 2018a; Mayourian et al., 2018b). Although our recent work has modeled paracrine signaling effects in these known cell types (Phillips et al., 2021), modeling the paracrine effects of custom cells was beyond the scope of this study and would first need to be determined empirically.

Third, we only provide single-cell metrics of arrhythmogenicity, as tissue-level analysis was beyond the scope of the current study. Herein, we show that custom hMSC and hCIC cells expressing cardiomyocyte inward rectifier channel had APD₉₀, APD₅₀, and therefore APD triangulation closer to healthy conditions. In addition, resting membrane potential better resembled healthy conditions, all favorable single-cell indicators of arrhythmias. Nevertheless, all coupled cell models have decreased upstroke velocity, which may indicate a pro-arrhythmic potential of custom cells. Arrhythmogenic risk of hMSCs has been addressed previously by our group using tissue-level analysis (Mayourian et al., 2016; Mayourian et al., 2018a) and a similar approach could be used with the accepted custom cells and/or engineered hMSCs and hCICs. Addressing this limitation may be clinically relevant, as the recent CONCERT-HF trial showed a higher incidence of ventricular arrhythmias in all cell therapy groups, with the highest rate of occurrence reported with isolated hCIC treatment (albeit these findings were not statistically significant) (Bolli et al., 2021). Notably, our simulations suggest hCICs expressing the cardiomyocyte inward rectifier channel led to increased upstroke velocity compared to control hCICs.

Fourth, we note that the stringent acceptance criteria could lead to false negatives; for example, it is feasible that the restoration of AP or CaT alone may be sufficient to provide meaningful physiologic benefit despite not meeting our predefined acceptance criteria for both AP and CaT. However, we chose to err on the side of being conservative in these exploratory studies. In addition, given the computational expense of implementing a genetic algorithm across many parameters, we imposed bounds on the searchable space for certain parameters (e.g., gap junctional conductance) based on previous findings (Valiunas et al., 2004; Phillips et al., 2021). However, the nature of this exploratory simulation work allows for supraphysiologic conditions, which could also lead to acceptable solutions. For example, accepted custom cell heterocellular gap junctional conductances were 5–8 nS; while this is higher than reported hMSC-cardiomyocyte coupling conductances (Valiunas et al., 2004), even larger conductances that are representative of cardiomyocyte-cardiomyocyte junctional conductances (Jongsma and Wilders, 2000; Dhillon et al., 2013) may also lead to acceptable solutions. Such parameters beyond our searchable space could be more favorable in terms of resulting tissue-level metrics such as conduction velocity, warranting further investigation.

Finally, we acknowledge that although mathematical models offer the capability to explicitly control variables that can be challenging to manipulate experimentally, they are also subject to the limitations of their underlying assumptions, and ultimately the most valuable insights will come from a synergistic combination of computational modeling and experimental validation. For example, herein we assume all cardiomyocytes will need to be sufficiently restored, requiring an abundant number of therapeutic cells to be delivered (2-5 custom cells per cardiomyocyte). However, recent gene therapy studies have shown that correcting only a small percentage of cardiomyocytes can improve muscle function (Long et al., 2016; Nelson et al., 2016; Dave et al., 2022). Further in vitro and in vivo validation studies are therefore warranted to explore the cell coupling efficiency required to achieve a therapeutically beneficial restoration effect.

To support the translation of our simulated work to the bench, we provide estimations of raw ion channel parameter values required to achieve restoration. There already exists a precedent of reprogramming non-excitable cells (Zhou et al., 1998; Chen et al., 2007), hMSCs (Bruzauskaite et al., 2016), and hCICs (Vegh et al., 2020) to express non-native ion channels that survive well, thus providing promise for a similar methodology. However, we acknowledge the practical challenge of titrating the expression of channels of interest to match the accepted range of solutions listed herein; for convenience, we provide practical reference data points for wet lab colleagues to validate our computational findings and potentially enhance cell-based cardiotherapies.