Making a model: Part 4 - manual changes (IKr)

Background

In the last section, we’ve seen how MGAs may be used to develop a model. However, when using them, it is quite possible you may get a plot of fitness of creatures in the population as follows (apologies for its ugly printscreen nature; it comes from a quick snap into a digital logbook, as is the case with most figures in this section):

On the x-axis is the fitness component corresponding to action potential morphology, on the y-axis is the compound value aggregating other criteria. You can see the tradeoff between the two dimensions easily. In case you’re wondering about what is the green rectangular object, it shows the zone that corresponds to what I thought would be a reasonably good solution at the time. Of course, nothing inside, nor even close :( This means one of two things; either there is a model parameter that we are missing in the creature DNA (good luck finding it among the hundreds of possible numbers in the model), or there is a structural problem in the model – an important current may be missing, or the formulation of a current is problematic.  Below, we’ll go over a concrete example that I encountered in model development of what one may do in the case where an automated optimizer fails to find a suitable solution. The abbreviation I selected for this process is CSOMSKY (Case Study Of Manual Searching, King in Yellow [1]).

I note that the plots used in this section are, unless specified otherwise, based on an intermediate development version of the model and are not representative of the final model version.

Diagnostics

The first step is to get proper diagnostics. The fitness function in Matlab returns only the fitness number[s], in our case two fitness dimensions. However, as mentioned previously, one of the dimensions was in fact aggregating more calibration criteria (calcium transient properties and response to I_Na+I_NaL, I_CaL, and I_K1 partial blocks). It’s thus useful to make a version of fitness function which also returns the calibration criteria in isolation, before they are weighted-summed into a fitness component. It may be also of use to return the traces of action potential, calcium transient, etc. in case one wants to inspect them too. This diagnostic-fitness function may be run after the MGA finishes, running it on creatures on the pareto front, gaining more information about them.

First steps

What I expected the most was to find a single calibration criterion that prevented achieving an overall satisfactory fitness. If that was the case, I’d have probably tried to find the difference between models which perform best and worst at the calibration criterion, taking it further from there.

Interestingly, the situation was different, and it was actually two criteria together causing the problem: one requiring the calcium transient amplitude to reduce in response to 50% sodium blockade (blocking both I_Na and I_NaL), and the other requiring APD shortening in response to 50% I_CaL block. Class A of creatures produced models which correctly reacted to sodium blockade, but incorrectly prolonged the action potential duration in response to the block of I_CaL. Class B of creatures was the opposite, manifesting an incorrect reaction to sodium block, but a correct one to the I_CaL block. There was a clear tradeoff between the two calibration criteria – some models were clearly good at one and horrible in the other, and then also models that were quite balanced, being slightly bad in both. For further investigation, I decided to focus on the models clearly in class A or in class B.

As the next step, I sought to find out what is the difference between the two classes; given that most of our creatures’ DNA was current multipliers, this can be worded as finding which currents differ mainly between the two populations [2]. The current most strongly and consistently different between creatures belonging to class A and B was I_NaL, the late sodium current (with class A having a high I_NaL).

Verification of the first step

It might still be very well possible that I_NaL is the main difference, but not the causative factor. To verify the causal role, I took creatures of class A (good sodium block, poor I_CaL block) and reduced their I_NaL. The creatures changed into class B, suddenly manifesting a poor reaction to sodium block, but a good one to I_CaL block. Thus, I_NaL appears to be a causative factor. For completeness of information, I measured the effect of I_NaL reduction in two ways which give slightly different information:

• Pre-pacing a control cell for 1000 beats, followed by acute I_NaL reduction for 5 beats and measurement of the cell properties during these. This gives information about the direct electrophysiological effect of the current reduction. In this particular case, the acute change induced a clear change in behavior, suggesting a direct electrical influence of I_NaL on the phenomena of interest.
• Chronic I_NaL reduction, i.e., reduction of I_NaL already from the start of simulation, followed by pacing for 1000 beats. This was then compared to a cell paced for 1000 beats without the I_NaL reduction. If there was no major change in cell behavior with acute change (which we did see, as described above), but there was an effect with a chronic change, it would have suggested the role of I_NaL is about a long-term effect such as a change in ionic concentrations.

Another approach that I did not use here, but that can be used to dissect causality, is clamping of variables. Do you suspect that CaMKII phosphorylation of a channel promotes EAD formation? Try to repeat the simulation with the said phosphorylation (on that site only) clamped not to exceed a given threshold – if that stops the EAD facilitation, it’s one more reason to think your suspicion holds water.

Finding the culprit

We have discovered that I_NaL plays a large role in the fork between good behavior of sodium and I_CaL blockade, however, this does not tell us anything about the mechanism. At the same time, finding the mechanism is critical, as that’s what can eventually lead us to an improvement.

I took a model of class A as a starting point, trying to understand why its reaction to I_CaL block is poor, showing an APD prolongation. Starting with class A model rather than class B was simply because a good reaction to sodium blockade was one of the most important criteria to us – it is in fact where the whole project on development of ToR-ORd started.

The AP morphology of the selected intermediate model in control condition (blue) and under I_CaL block (red) looked as follows:

(x axis is time - last spike of many, hence the weird range, and y axis is the membrane potential)

For this part of investigation, I used a 100 % I_CaL block to make the effect as clear as possible. Subsequently, I plotted the currents in control and I_CaL block conditions, looking for main differences in these.

Unsurprisingly, the largest difference in currents was in the reduction of I_CaL with its full block, but this would strongly shorten the APD if other factors were kept constant, so it was certainly not the explanation. Rather, it was clear that there are other factors which override this trend. The two currents changed most prominently by the full I_CaL block were increased I_NaL, and reduced I_Kr, both contributing to the APD prolongation.

Culprit 1 – I_NaL

The I_NaL in the two conditions (identical color coding to the previous plot) looked as follows:

The increase in I_NaL amplitude following I_CaL block could be due to changes in activation, inactivation, or driving force. An inspection of state variables during an action potential revealed that the dominant factor in the I_NaL increase is the increase in driving force, which follows from lowered plateau potential under I_CaL block.

Is such an increase in I_NaL plausible? In general, yes, the increase in driving force is hard to argue against. However, one may wonder if there aren’t other factors coming into play, linked to special gating properties of the Nav 1.5 channel mediating I_Na and I_NaL. Studies based on single-pulse voltage clamp predicted that I_NaL would be shaped like an extended triangle-like decay of I_Na; however, a recent-ish study (Horvath et al. 2013) employing an AP clamp technique suggests it is not really the case (Figure 1A in this article), explaining the difference via non-equilibrium gating (there are some very interesting papers on this, e.g. (Clancy et al. 2003))[3]. If it is the case that I_NaL is generally small shortly after AP upstroke as shown in the Horvath et al. paper, perhaps the effect of increased driving force is mitigated. However, in the absence of data to decide what is happening with I_NaL under I_CaL block, it was hard to say what to do and I decided to switch focus to the considerably reduced I_Kr.

Nevertheless, at this stage, it dawned on me why I_NaL was at the core of the fork between classes A and B in the population returned by the MGA. As described in the ToR-ORd paper, a reduction in I_NaL reduces the calcium transient amplitude; mainly via APD shortening and resulting smaller calcium influx via If I_CaL. If I_NaL is generally high, its reduction causes a more prominent APD shortening and thus a greater reduction in calcium transient amplitude. This is why models with higher I_NaL could manifest a net calcium transient amplitude reduction in response to 50/50 I_Na+I_NaL block. However, at the same time, given that I_NaL is increased under I_CaL block due to the increased driving force, models with high I_NaL show a greater total increase in I_NaL following an I_CaL block, and thus manifest APD prolongation, performing poorly on the criterion of I_CaL block and APD shortening. Despite the utility of understanding the fork between classes A and B, it was pretty clear that yet again, understanding a problem does not mean it’s obvious how to solve it.

Culprit 2 – I_Kr

The second factor involved in the APD prolongation in response to I_CaL block was I_Kr; at the time, the current formulation was the same as in the O’Hara-Rudy (ORd) model. Following  I_CaL blockade, I_Kr was reduced considerably:

Again, this reduction could be, in theory, due to changes in inactivation, activation, or driving force. Inspection of state variables has shown it is the activation – and it could be either a change in the steady-state value of activation, or the time constant of reaching it. Plotting the steady-state values for the given membrane potentials versus the actually achieved membrane potential showed a marked difference, pointing to the time constant as the causative factor.

The reason for the lack of I_Kr is described in our paper, but to briefly recapitulate it here, let’s look at the plot of the fast and slow activation time constants of ORd model’s I_Kr (to the left, taken from the original ORd article), and the action potential morphology of the final ToR-ORd [4] as well as ORd (to the right, taken from our paper):

The membrane potential during the early plateau in our model (both the intermediate and final ToR-ORd) is around 20-23 mV, and this clearly corresponds to low values of fast and slow I_Kr activation time constant, and thus relatively rapid activation. However, when the early plateau is reduced by I_CaL blockade down to ca. 0 mV, we can see that the time constants become very high (more than a second), meaning that the I_Kr does not have enough time to activate to the steady-state value. At the same time, the same problem does not manifest in the original ORd model – given its very high plateau potential, a reduction  in plateau caused by I_CaL blockade does not increase I_Kr activation time constants drastically.

Having understood the problem, what do we do about it? [5]. Going over various publications on I_Kr activation, I noticed that data underlying ORd I_Kr have generally fairly slow activation dynamics at the membrane potentials of -10 to 20 mV. Trying to reduce the activation time constants ameliorated the problem somewhat, but was not enough. Simply shifting the activation time constant curves to the left would be another approach we considered, but to see a sufficient difference, the shift would have to be 15-20 mV, which is a lot. It would be quite debatable [6] if the model might be called data-driven at that stage without further justification. I then started exploring more complex Markov models of I_Kr, which looked vaguely promising. A complication to the whole process of I_Kr explorations was that we had no direct evidence that the ORd I_Kr model is problematic - it could be that the rest of the model needed fixing.

The breakthrough came when my DPhil supervisor Prof. Gil Bub was visiting Oxford, and with him his PhD supervisor, Prof. Alvin Shrier. Going for a pint of beer together in the evening, I explained the problem and Prof. Shrier suggested I should look at the work of Prof. Jamie Vandenberg on I_Kr and explained how the single-pulse voltage clamp protocols often fail to characterize I_Kr dynamics well, e.g., mixing activation and recovery from inactivation together.

A question that might initially look like a sidestep – what might one hope for in the wildest optimistic dreams? I suggest the following:

• The Vandenberg group published a model of I_Kr that does not manifest the major reduction in current with reduced plateau following an I_CaL blockade.
• The model is easy to incorporate and works without major changes either in itself, or in the rest of the model.
• The paper would also conveniently publish data demonstrating that lowering the AP plateau should not produce a drastic reduction in I_Kr.

Guess what – the paper on the Lu-Vandenberg model of I_Kr does just that! Moments like the one when I plugged the Lu-Vandenberg model in and it worked out of the box are the reason why science is so addictive. Looking back at the meeting in the pub, I guess that corroborates the ages-old wisdom In vino veritas.

Take-home messages

Key points in my view are:

• The fact that a model is data-driven does not mean that it is bulletproof or that it should be taken as set in stone. For example, as seen in literature on late sodium current and I_Kr, conventional single-pulse voltage clamp protocol may not produce sufficient information to characterize the behavior of a current when the gating is complex.
• Linking a model’s behavior to a change in a particular current density is a good start of investigation, but seldom the answer. What is described in this blog post started with a fork between good/bad behavior under block of sodium and calcium currents respectively… and the answer was in the replacement of a potassium current formulation. I personally would never discount the possibility that that a computer model gives us answers we like [7] via a problematic mechanism, and so it is important to get to the bottom of the mechanism whenever possible, rather than just associating behaviours with currents.
• Understanding mechanisms is also invaluable in model development whenever an obstacle is reached – it doesn’t give an easy how-to on how to improve the model, but without it, it’s even harder.

---

[1] The King in Yellow is present only to make the abbreviation work nicely – nevertheless, it is an interesting and not that well known old book. Especially if you’re fans of E. A. Poe or H. P. Lovecraft, give it a try – it’s freely and legally available via Project Gutenberg.

[2] It’s good to keep in mind the total current amplitude/importance. If there is a relative change of 20% in a major current like I_CaL, it may be a more important difference than a 50% change in the background calcium current.

[3] That said, the AP clamp is also not fully bulletproof approach, as the current blocker used may not be fully specific (it may also block yet unknown currents/processes), and the current blockade may elicit changes in ionic concentrations that also modify the cell behavior. I’d say one rather clear case is when L-type calcium current is estimated as nisoldipine-sensitive current - it is reasonable to expect that blockade of calcium influx will have a major effect on the current carried by sodium-calcium exchanger, and this change will be also reflected in what is measured as the nisoldipine-sensitive current.

[4] We worked with an intermediate development version of our final model in this section so far, but for the illustration of the point here, the final ToR-ORd AP shape can be used, given its similar plateau potentials to the intermediate model.

[5] This may sound like a common theme of the model development at this stage. However, this is still a relatively optimistic phase. It’s reasonable to expect spending more time in the state where one does not understand what the problem is in the first place, let alone having an idea of the solution.

[6] When reading this sentence, please take into account that I live in the UK for more than 5 years. For more direct styles of communication, the translation to international English would be more like “No, this really is not data-driven modelling.” It took me a while to understand that this funny table is an accurate depiction of reality and when I stopped taking it as a joke, I found communication with native English persons much easier to navigate.

[7] I suppose this is linked to a well-known but hard-to-guard-against bias, where one very critically approaches a problem when things don’t work (looking for new ways to analyze data, model, etc.), but when the data analysis fits the expectations, the rigor of criticism suddenly drops and a paper starts being written. But precisely because it’s a very natural thing and it’s convenient to succumb to this bias, it’s so important to be aware of it and combat it vigorously. As also seen in online debates on anything more serious than kittens, I think this is not about education/knowledge – it’s about attitude and self-criticism.