Representation of Multiple Cellular Phenotypes Within Tissue-Level Simulations of Cardiac Electrophysiology

Abstract

Distinct electrophysiological phenotypes are exhibited by biological cells that have differentiated into particular cell types. The usual approach when simulating the cardiac electrophysiology of tissue that includes different cell types is to model the different cell types as occupying spatially distinct yet coupled regions. Instead, we model the electrophysiology of well-mixed cells by using homogenisation to derive an extension to the commonly used monodomain or bidomain equations. These new equations permit spatial variations in the distribution of the different subtypes of cells and will reduce the computational demands of solving the governing equations. We validate the homogenisation computationally, and then use the new model to explain some experimental observations from stem cell-derived cardiomyocyte monolayers.

Keywords: Bidomain; Cardiac electrophysiology; Homogenisation; Monodomain; Stem cell-derived cardiomyocytes.

Fig. 1

Simulated action potentials of atrial-like and ventricular-like human stem cell-derived cardiomyocytes, generated using the Paci et al. (2013) model. Two properties of the action potentials are indicated. The maximum diastolic potential, MDP, is the most hyperpolarised potential. The action potential duration, ${\text{APD}}_{90}$, is the time taken to achieve a given percentage (here, 90%) of repolarisation following the upstroke

Fig. 2

Different spatial distributions of cellular phenotype. The case on the left may be divided into two partitioned regions, each containing a single type of cell. Partitioning the case on the right into single-phenotype regions would result in many tiny partitions. Performing the homogenisation process over regions containing both types of cell is therefore preferable

Fig. 3

Spatial layout of the different phenotypes in the first three sets of simulations. The dark and light shades represent regions within which one of two cellular phenotypes is exclusively present. The cellular electrophysiology models that represent the two cell types are chosen from the six parameterisations of the FitzHugh–Nagumo model described in the main text. Intermediate shades denote the HP model with appropriate values of ${\rho }_1$ and ${\rho }_2$ (the relative contributions of each phenotype). The value of n indicates the number of regions into which the domain was partitioned when the PP model was used

Fig. 4

Spatial layout of the different phenotypes in the final three sets of simulations. As in the previous figures, the dark and light shades represent two different phenotypes. In Sets 4 and 5, we utilise the parameterisations of the FitzHugh–Nagumo model that are listed in the main text, while in Set 6 we use atrial-like and ventricular-like models of hSC-CM electrophysiology (Paci et al. 2013)

Fig. 5

Action potentials of the six parameterisations of the FitzHugh–Nagumo model. The three self-exciting models (top) beat at their natural frequencies, while the excitable models (bottom) are stimulated every 500 time units. Upstroke times have been aligned at $\text{time}=0$

Fig. 6

Variation in ${\text{APD}}_{90}$ (time to achieve $90\%$ repolarisation) during the final beat in Set 1 simulations. The panels on the left show the variation in ${\text{APD}}_{90}$ across the fibre during the final complete beat for three selected cases: those of the PP model with the largest and smallest partitioned regions, and the HP model. The panels on the right show the minimum and maximum values of ${\text{APD}}_{90}$ across the central region of $35<x<65$ during the final complete beat of all Set 1 simulations. Values from the PP model are shown using crosses, while those from the HP model are indicated with the dotted line

Fig. 7

Variation in maximum diastolic potential (MDP) from Set 1 simulations. As in the previous figure, the panels on the left show the MDP across the entire fibre for selected cases. Panels on the right provide a summary of the minimum and maximum values of MDP recorded in all simulations

Fig. 8

Conduction velocity of the travelling waves in the Set 1 simulations. Activation times in cases where the HP model was used to simulate a spontaneously activating system (i.e. models S1–S3 and S1–E2 in the top panels) were synchronous, leading to an infinite conduction velocity

Fig. 9

Beat rate from Set 2 simulations. The beat rates of fibres simulated using the PP model are compared to those of fibres simulated using the HP model with equivalent proportions of the two phenotypes. With both sizes of partitioned units, the discrepancies between the HP and PP models (indicated by the proximity of the cross-dot pairs) are generally small. The discrepancies are noticeably smaller in the lower panel, where the smaller partitioned units are used. The major difference between the small and large partitioned unit simulations may be seen around ${\alpha }_H=0$. The HP model is quiescent at this value of ${\alpha }_H$, as are all instances of the PP model with small partitioned units. However, spontaneous beating is still seen in some of the simulations that utilise the PP model with larger partitioned units (Color Figure Online)

Fig. 10

Activation times of all recorded beats in PP model simulations (Set 3) with equally sized partitions. The activation times are normalised so that the earliest activation time during each beat is set to 0. Model S1 has a slower natural beat rate than Model S3, and so propagation spreads from regions where Model S3 can dominate. As the number of partitions increases, the apparent conduction velocity increases due to synchronisation effects as we tend to the homogenised case (Color Figure Online)

Fig. 11

Activation time of the final beat for Set 4 simulations. We compare the homogenised phenotypes (HP) model (black line) to the partitioned phenotypes (PP) model with 60, 120 and 240 partitioned units. Parameter A controls the extent of variation in phenotype—with more variation in phenotype there is a slower wave speed in both HP and PP models and a more noticeable difference between the HP and large-unit PP models. See Fig. 4 for the underlying phenotype arrangements across this domain

Fig. 12

Activation times of the final beat in the Set 5 simulations. The bold black line shows the activation time from the homogenised phenotypes model. The thinner grey lines show the activation time from the 15 randomly assigned partitions for the PP model [probabilities given by Eq. (12)], so that each has a slightly different phenotype layout, examples of which may be seen in Fig. 4. The HP model and PP model wave speeds are in good agreement (the same gradients are seen in these activation time plots), but the random arrangement of phenotype partitions can change the location of the emergent ‘pacemaker’ site(s) in the PP model

Fig. 13

Activation times of the final beat of simulations in the Set 6 simulations. The two cell types are represented by the Paci et al. (2013) models of ventricular-like and atrial-like hSC-CM electrophysiology. Results from the homogenised phenotypes model are shown as a bold black line; the thin grey lines represent the 15 randomly generated partitions according to Eq. (12). Each repeat has a slightly different phenotype layout; see Fig. 4 for examples of the underlying phenotype arrangements across this domain. As we noted in the previous set of simulations, the HP model and PP model wave speeds are in good agreement and we observe that the random arrangement of phenotype partitions can change the emergent ‘pacemaker site’ in the PP model

Fig. 14

Cells of two different phenotypes in the HP model. The repeated unit is made up of both cell types; it should be noted that the two types do not have to be present in the same proportions. ${\Omega }_i$ and ${\Omega }_e$ denote the intracellular and extracellular domains, respectively. The surface area of the cell membrane within the repeated unit is given by ${\Gamma }_m={\Gamma }_{m1}+{\Gamma }_{m2}$. The lengthscale of the solution, L, is assumed to be much larger than the lengthscale of the repeated unit, l