Making a model: Part 5 - manual changes (ICaL)

A second topic of biggish advancement is the changes to L-type current (I_CaL) pertaining to the ionic activity coefficients and activation curve (all described in ToR-ORd paper), which this blog post is about. It’s much less of a step-by-step compared to the previous post, as it was much more about conditioning the brain and waiting for a heureka moment.

The spark that triggered ToR-ORd development was when I was playing with some published data on post-infarction remodelling in ORd (including a major reduction in the fast sodium current I_Na) and got a major increase in calcium transient amplitude instead of the expected reduction. I tracked it down to the reduction in fast sodium current among all the remodelling factors, which alone was causing a major increase in the calcium transient amplitude. The next step was to go through the literature on sodium current and its pharmacological reduction, which suggested rather unequivocally that reduction in sodium current is supposed to reduce the calcium transient amplitude.

Going over what happened in the ORd when I_Na was reduced, it was clear that the increase in calcium transient was due to a major increase in I_CaL[1]. This was particularly pronounced around the peak of I_CaL, and shortly after, and given than in the model, SR release directly depends on I_CaL, this translated into an increased SR release and thus calcium transient. This increase in I_CaL could be in theory due to 1) driving force, 2) activation, 3) voltage inactivation, 4) calcium inactivation.

Starting with inactivation, it was fairly clear that reduced degree of I_CaL voltage inactivation upon I_Na block was contributing to the I_CaL increase, but this is logical – with reduced peak and very early plateau following I_Na block, the reduction in voltage inactivation is natural. I.e., this was unlikely to offer much room for remedy of the discrepancy with experimental data. I also did not think it likely that calcium inactivation would be of much help – a stronger calcium-driven inactivation could limit I_CaL increase perhaps (because increased I_CaL leads to more SR release and thus stronger calcium inactivation), but this would require increased SR release in the first place, which is something I was trying to get rid of.

Next I focused on the driving force, learning about the Goldman-Hodgkin-Katz (GHK) flux equation used for L-type current since the first Luo-Rudy model of 1994. There, it was critical to have the opportunity to discuss this with Prof. Yoram Rudy himself who was visiting Oxford at the time. He sent me extremely useful literature on the Goldman-Hodgkin Katz equation and the Debye-Hückel theory of driving force on which the driving force of I_CaL since the Luo-Rudy model was based.

The driving force based on the GHK equation is:

where z is the charge of the given ion, V is the membrane potential, F,R,T are conventional thermodynamic constants, and [S]i, [S]o are intracellular and extracellular activities of the given ionic specie. The activity S is computed as a product of the activity coefficient and the concentration (in either the intracellular or extracellular space). The activity coefficients are usually based on the Debye-Hückel equation (See Wikipedia for an introduction) or one of its variants (Guntelberg or Davies equation).

In general, the equations made good sense to me, except I struggled to find how the ionic activity coefficients were derived (0.341 for extracellular space, 1 for intracellular). Ionic activity coefficients are used to “weight” the intra/extracellular concentrations in the equation, where the fewer ions are around in a solution, the higher the activity. It can be illustrated, e.g., on humans fleeing a house upon learning the house is full of vampires, where the house has relatively narrow exit doors (people corresponding to ions and doors to a channel). If there are only few people inside, they can move freely and they can pass through the exit with ease. However, if the house is crowded, people start running and bumping into one another, or even fighting to get to the exit (corresponding to interaction between charged particles). Consequently, the rate of people leaving the building (~ionic flux) is lower than one might expect from so many people wanting to leave the house.

However, no matter whether I used the Debye-Hückel equation itself or one of its extensions, I could never get the values of 1 and 0.341. The next month or so was spent chasing various colleagues who even barely touched chemistry or physics (not one of them has heard of the GHK equation). Then I even went to library to go through ancient tomes of wisdom. After several days of this, I finally found a book about thermodynamics of membranes, where the GHK equation was presented, yay! Or, rather, nay, because when I opened the corresponding pages, only the basic version of the GHK equation was present, assuming the ionic activity coefficients to be both 1.

Feeling quite unhappy about failing to resolve the mystery quickly, I realised at some stage how strange the coefficient of 1 for intracellular space is. Based on any of the Debye-Hückel equation variants, this is possible only when there are no ions. There I saw it might be not my head where the problem is, but that the coefficients themselves might not be accurately estimated. Picking various bits of information from the library books and from online resources, it seemed like the intracellular and extracellular ionic activity should be pretty similar, particularly because the ionic strength is likely to be similar. I thus chose Davies equation, a variant of Debye-Hückel equation suitable for the relatively high ionic strength of solution inside and outside myocytes and used that, getting ionic activity coefficients of between 0.6-0.7 for intracellular and extracellular space [edit - see an update on this here]. I also decided to compute the coefficients dynamically at each point of the simulation as to reflect changes in ionic concentrations – there isn’t a major simulation time penalty for this. I still wonder about how the original values were obtained, but the current hypotheses I have are not safe for mortal’s eye and might be best kept in a Supplement to Necronomicon until the day when the stars of the Old ones enter a conjunction (basically, I don’t really have a good idea). Another intriguing historical mystery is that the Shannon et al. rabbit model of 2004 claims to have used the same approach as the LRd model with regards to driving force, however, it used 0.341 for both intracellular and extracellular space. This has subsequently propagated into other popular models, such as by Grandi et al. 2010 Was it an accident or intentional change? In a way, this is, I’d say, a more suitable combination of ionic activity coefficients than 0.341&1. The absolute values are not particularly important (the driving force is then scaled by I_CaL conductance anyway, so whether both coefficients are around 0.34 or around 0.6 doesn’t matter), but whether the intracellular/extracellular activity coefficients are similar or not does play a role (e.g., in current reversal). Going back to the original question of I_Na block and I_CaL increase, the change in activity coefficients seemed to have helped a bit, but not quite enough…

A second important factor in the I_CaL update (one that is also relevant for most existing models, I believe) was the realisation that not all is well with the activation curve [2]. First, I was thinking in the direction “this is the last thing we can change, and if we make it steeper between peak potential of control and I_Na-blocked model, maybe we can ameliorate the I_CaL hyperactivation upon I_Na block”. To this end, I searched the literature for the steepest activation curve in that interval, finding one (Harasztosi et al. 1999) in a debatably relevant type of tissue (myotubes derived from skeletal muscle). This seemed to again ameliorate the issue at least partially. However, I wasn’t overly happy as the dataset was an outlier in this regard, with more human cardiomyocyte datasets looking more like the Magyar et al. 2000 used in the ORd model. I appreciate there is a lot of biological and measurement-induced heterogeneity, but if the only way of achieving a decent model is to rely on outlier data from only vaguely relevant cell type, I don’t think that’s ideal and massively convincing.

The problem sorted itself fortunately, with the realisation that if we use GHK equation to compute the I_CaL based on its activation (among other things), the activation curve should be computed as the experimentally measured IV curve normalised by the GHK-equation-based driving force. I then went into the literature to find which driving force was used in the articles from that era, acquiring a strong impression that it was simply based on the V-E_Ca equation with experimentally measured reversal potential E_Ca of 60 mV, rather than the GHK equation. Manually extracting the points on the I-V relationship and activation curve in the work by Magyar at al., trying to reconstruct the driving force used to get the activation, it became crystal clear that it was indeed V-60 used as the driving force. Again, I wasn’t sure for a while whether this inconsistency was my misunderstanding of literature, or if it was real (with all models relying on the GHK equation not having this right). What sort of decided it for me was the fact that when I replaced the original ORd activation curve with one based on I-V relationship normalised by an estimate of the GHK, the ok-but-not-great result of simulation of I-V relationship in the model suddenly became pretty much spot on (Figure 2D in the ToR-ORd paper), which was like an internal mini-validation. I then presented it to my PI (Prof. Blanca Rodriguez) and the group and it seemed that everyone agreed that this was the way to go. Beyond the improved I-V relationship, one immediately appealing aspect of the new activation curve was its steepness in the zone of membrane potentials affected by I_Na block, and after this change, I_Na block finally lead to only a very meagre increase in I_CaL – I finally didn’t have to resort to outlying myotube!

In parallel with the developments of driving force and activation, it occurred to me that “sodium blocker” drugs do not block only I_Na, but also I_NaL (many even preferentially). This is something I retrospectively kick myself for, as I should have realized it sooner, but I guess one cannot have it all. I_NaL prolongs APD and thus extends the duration of I_CaL, enhancing total influx of calcium into the cell, which increases the sarcoplasmic reticulum calcium loading. Consequently, increased I_NaL increases calcium transient. When I tested an unbiased 50/50 block of I_Na and I_NaL, at long last, the net effect was a reduction in calcium transient amplitude – case closed (of course only for a day or so, before I found out about the problem described in the previous section).

All these changes reshaped the I_CaL slightly, mainly shortening it given a different shape of activation curve. Consequently, the changes shortened the APD, and in such a way that reducing e.g., I_Kr to counterbalance the shortening did work, but the resulting AP shape wasn’t great. To combat this, all that was needed was to recall literature showing that a part of I_CaL channels is present outside the dyadic subspace. I assumed this “main membrane I_CaL” would last longer than the subspace one given less calcium-dependent inactivation, and such was the case – upon placing 20% of I_CaL outside the junctional subspace, the AP morphology started looking nice and consistent with data again. It's an interesting semi-philosophical question of how appropriate it is to make such choices based on needing to get a nicely-behaved cellular model. It can be easily argued that we should include such knowledge if we know of it (this doesn't even make the model more complex), whether or not it helps us. 

Retrospectively, all these observations/insights are quite straightforward [3], but they weren’t to me before they materialised, probably given the many alternatives that had to be considered. Also, there was no really clear path to realizing the insights – most of them happened suddenly when I was relaxing with tea or in nature, hovering at the edges of consciousness, waiting for bits of knowledge that were chaotically flying in my head to crystallize. Anyway, the take-home message here is: Don’t get into the frame of mind “I cannot be right with this interesting observation, because if it was true, someone else would have found it by now.” I have heard people say this and it is an unnecessary self-demoralising limitation. Yes, it’s definitely good to err on the side of caution and be modest, but if you cannot invalidate your observation with reasonable effort, it’s worth accepting you may be right and it’s at least worth bringing the observation to someone else for consultation.

---

[1] The problem was much less pronounced in the canine Heijman-Rudy et al. (2011) model, which is why the model was suitable for a previous studies on infarct remodeling on which I worked (Tomek et al. 2019, Tomek et al. 2017). It was quite a lucky choice, as I was unaware of the whole issue of sodium currents and contractility at the time. I tried to port the canine I_CaL model to ORd, changing its properties (activation curve etc.) to fit the human data, but I could never make it work all-around, either not being able to get EADs, or not getting the reduction in calcium transient amplitude upon sodium blockade. I think among the key reasons why the Heijman-Rudy didn’t manifest an excessive increase in I_CaL in response to simulated sodium blockade is its activation curve. While ORd’s original activation curve was basically flat from 15 mV on, Heijman-Rudy I_CaL model flattens its activation only around 50-60 mV. Given that I_Na block in Heijman-Rudy model considerably reduces peak of action potential, this reduces I_CaL activation, and the peak current then does not increase in the fashion that ORd’s does. However, the activation curve of Heijman-Rudy I_CaL model is different from human data (which is not a criticism, by the way – it’s based on animal data), and when I replaced that with the human-derived one, the model stopped reacting to sodium blockers well. Another important factor is that the sodium blockade in ORd versus Heijman-Rudy (or ToR-ORd for that matter) is the AP morphology. In ORd the I_Na block affects also early very plateau potentials (see Figure 3B of our ToR-ORd paper), which makes the effect of the block on I_CaL driving force and voltage inactivation last longer, compared to Heijman-Rudy or ToR-ORd, where the early plateau is relatively unchanged and the main difference is around the action potential peak, which lasts a shorter time.

[2] This is obtained as the experimentally measured I-V relationship (how much I_CaL current flows in response to pulses of increasing potential) divided by the driving force at the given membrane potential.

[3] A quotation comes to my mind, this time from our professor of mathematical structures and algebras, Prof. Aleš Pultr: „There are only two types of problems: the ones that are trivial, and the ones that are yet unsolved.“