Building a model from a publication

All the equations are there, you’ve got figures of model voltage traces that show what the model is supposed to look like, and you’ve got a table of parameters that the authors used to get those figures. It might seem like a breeze to type all that up in Matlab and pop out identical results in reproducing that conductance-based model, but your patience and methodicalness will be rewarded. I recommend starting with the function code that will be called by your ode solver and building in one conductance at a time (see example below).

function dy = myNeuron(t,y)
Vm = y(1); % mV
H = y(2); % gating variable for a conductance
M = y(3);

cm = 1; % nF
g_leak = 0.001; % microS
g_k = 0.2;
g_Na = 0.3;
v_leak = -40; % mV
v_k = -80;
v_Na = 100;

h_inf = 1/(1 + exp((Vm + 15)/6)); % hypothetical steady state gating
m_inf = 0.5*(1 + tanh((Vm - 31)/22));

I_leak = g_leak*(Vm - v_leak); %nA
I_K = g_k*H*(Vm - v_k);
I_Na = g_Na*M*(Vm - v_Na);

dy(1) = (-I_leak - I_Na - I_K)/cm;
dy(2) = (h_inf - H)/tau_h;
dy(3) = (m_inf - M)/tau_m;

This could be your inner function, myNeuron, that is called by the ode solver you choose.

[T,Y] = ode45(@(t,y) myNeuron(t,y),[0 1],y0);

Start with the leak conductance alone – making the most boring, passive neuron ever.

Make sure that no matter what your initial voltage is, the voltage trace that you simulate will relax to the resting voltage you’ve chosen. Does it look like it’s taking too long to get there? Or getting there too quickly? Now is a great time to check your units and the time scale or time step you’re using. If anything about this first step of simulating a passive neuron got weird:

• Make sure the units check out. Sometimes authors publish their parameters in terms of conventional units but the units wouldn’t check out on both sides of the equation. For example, without putting a 1e-3 in front of the nA. Below is my little cheatsheet of unit conversions.
• Check the implementation of the ode solver you’re using. Make sure the variables haven’t gotten switched, and make sure they’re being assigned to the solver’s channels. Matlab documentation is helpful here: Example of Implementing ODE Solver and Parametrizing Matlab Functions.

unit	conversion
Voltage (Volt)	Current * Time / Capacitance (Ampere * second / Farad)
Current (Ampere)	Charge / Time (Coulomb / second)
Capacitance (Farad)	Charge / Voltage (Coulomb / Volt)
Conductance (Siemen)	1 / Resistance = Current / Voltage (1 / Ohm) = (Ampere / Volt)

giga	1e9
mega	1e6
kilo	1e3
unit	1e0
centi	1e-2
milli	1e-3
micro	1e-6
nano	1e-9
pico	1e-12

If you’ve got a working passive neuron with only a leak conductance, you can start building in the rest of your conductances. I recommend dividing your conductances into “major” and “minor” conductances. The major currents are leak, one main inward, and one main outward current. Work with those first to make sure you get something reasonable that qualitatively resembles a less complex model with only “major” conductances. For example, the Hodgkin-Huxley model only has leak, sodium (Na+), and potassium (K) conductances; or the Morris-Lecar model which has only leak, calcium (Ca), and potassium (K) conductances. Then add in your other currents one at a time and observe that they are changing the model behavior in a way that is expected. For example, including another K conductance should slow down spiking or eliminate it while adding in a Ca conductance should make a more excited response or produce some slow-wave bursting dynamics.

Troubleshooting a complete model

If you’ve gone through the above steps, or if you came to this guide with a complete model that isn’t working. Here are my tips for fixing it.

• Stupid things like typos are common. It might not even be your typo but the author’s typo. If possible, check the paper you’re using against another paper that uses the same model. Look for any discrepancies. Generally, the conductance equations should be exactly the same. It’s rare that those change from paper to paper. The parameters would be different from paper to paper (or even within the same paper), but make sure they look reasonable and use comparable units (i.e. it would be weird if one paper used 0.1 picoS for gNa but another paper used 100 mS since the difference is several orders of magnitude and it’s unlikely that both models are functional. Usually one of those is a mistake).
• The ode solver may not be right for the model. If you have lots of dynamical variables (>3), the likelihood of stiffness might go up since that can happen with variables that change at very different time scales. One way around this is to choose a very tiny time step while troubleshooting. A tiny time step will eliminate the problem of stiffness for any ode solver you use, so even though it will take much longer to simulate a short run, you can rule out other stuff in the meantime. 
• Or the simulation time scale is too coarse in general and you’d need to simulate on a finer timescale anyway even if you’ve chosen the best ode solver for your system.
• Go through each conductance and do a sign check. Make sure that the conductance produces a current of the correct sign and has qualitatively appropriate dynamics (fast, slow, etc.). Do this with the other currents off. Make sure the current is changing the voltage in a way that seems reasonable/expected (i.e. an inward current should depolarize your voltage). 
• Look for interdependencies in the conductance equations. If there are any, test those conductances together and make sure one isn’t causing the other to be dysfunctional.

If all else fails, I’ve had success with simply starting from scratch with clean code (no copy/pasting sections of the old code). I know that seems daunting, but you’ll wish you had done it sooner if it solves your problem. I’ve done this before and never ever found the mistake in my original code even though the new code actually worked – so sometimes a small error is really hard to find even when you know it’s in there somewhere.