We want to explore here the basic properties of the extra-cellular potential generated by a uniform active—that is, able to propagate an action potential (or spike)—cable or axon. We are going to use the conductance model of Hodgkin and Huxley (1952) together with the cable model making the "full" H & H model.
Remember that H & H did not solve their full model in their opus magnum, remember also that the mechanical calculator they had at this time was far less powerful that any of the smartphones everyone has nowadays in his/her pocket. They used a "trick" looking at the propagation of a waveform without deformation at a constant speed (\theta), that is, a spike or action potential. In this way the spatial derivatives of the membrane potential can be expressed as time derivatives (as we will see bellow) and the partial differential equation (PDE) of the full model can be replaced by ordinary differential equations (ODE).
So we are going to start by deriving the expression of the extra-cellular potential generated by a "cable like neurite"—a neurite with a large length to radius ratio—that can approximated by a line source. We will follow a classical development that is very clearly explained in the book of Plonsey and Barr (2007) Bioelectricity. A Quantitative Approach, published by Springer. This development will lead us to an equation relating the extra-cellular potential to the integral of the weighted second partial derivative of the membrane potential with respect to space, (\partial^2 V_m(x,t) / \partial x^2). The H & H model will give us actual values for this derivative but that will require a numerical solution. We will then explore the effect of the axonal diameter on the extra-cellular potential.
The electrostatic potential (\Phi_e) [mV] generated by a constant point source of intensity (I_0) [mA] is given by:
\begin{align}\label{eq:stat}\tag{1} \Phi_e = \frac{1}{4 \pi \sigma_e} \frac{I_0}{r} , ,\end{align}
where (\sigma_e) [S/cm] is the conductivity of the extracellular medium assumed homogeneous and (r) [cm] is the distance between the source and the electrode (Plonsey and Barr, 2007, Bioelectricity: A Quantitative Approach, p. 29).
For an extended source with a large length to diameter ratio (a cable) that can be approximated by a line source; the generalization of the previous equation for a continuous line source along the x axis (between (x_{min}) and (x_{max})) of the 3D Euclidean space equipped with Cartesian coordinates when the electrode is located at ((X,Y,Z)) is:
\begin{align}\label{eq:stat1}\tag{2} \Phi_e(X,Y,Z) = \frac{1}{4 \pi \sigma_e} \int_{x_{min}}^{x_{max}} \frac{i_m(x)}{r(x)} dx , ,\end{align}
where:
\begin{align}\tag{3} r(x) \doteq \sqrt{(x-X)^2+Y^2+Z^2};,\end{align}
and (i_m(x)) [mA/cm] is the current density at position (x) [cm] along the cable.
We get an expression for (i_m(x)) by considering a small piece of cable of radius (a) [cm] and of length (\Delta x) [cm] (Plonsey and Barr, 2007).
If the intracellular potential at position (x) is written (\Phi_i(x)), then Ohm's law—the current equals the potential drop multiplied by the conductance—implies that the axial current (I_i(x)) [mA] is given by ((\sigma_i) [S/cm] is the intracellular conductivity):
\begin{align} I_i(x) &= -\pi a^2 \sigma_i \frac{\Phi_i(x+\Delta x) - \Phi_i(x)}{\Delta x} \nonumber \ &\xrightarrow[\Delta x \to 0]{ } -\pi a^2 \sigma_i \frac{d \Phi_i(x)}{dx} , . \label{eq:stat2}\tag{4} \end{align}
Then the charge conservation implies that the membrane current density (i_m(x)) (positive for an outgoing current) is given by:
\begin{align} I_i(x+\Delta x) - I_i(x) &= -i_m(x), \Delta{}x \nonumber \ \frac{d I_i(x)}{dx} &= -i_m(x). \label{eq:stat3}\tag{5} \end{align}
Combining equation 4 and equation 5 we get:
\begin{align} \label{eq:stat4}\tag{6} i_m(x) &= \pi a^2 \sigma_i \frac{d^2 \Phi_i(x)}{d x^2}, . \end{align}
Now, writing the membrane potential (V_m = \Phi_i - \Phi_e) we have:
\begin{align} \label{eq:stat5}\tag{7} i_m(x) &= \pi a^2 \sigma_i \frac{d^2 V_m(x)}{dx^2} ,. \end{align}
This allows us to rewrite equation 2 as:
\begin{align}
\label{eq:stat6}\tag{8}
\Phi_e(X,Y,Z) = \frac{a^2 \sigma_i}{4 \sigma_e} \int_{x_{min}}^{x_{max}} \frac{1}{\sqrt{(x-X)^2+Y^2+Z^2}}
\frac{d^2 V_m(x)}{dx^2} dx ,.
\end{align}
The quasi-static approximation (Plonsey, 1967, The bulletin of mathematical biophysics 29:657-664; Nicholson and Freeman, 1975, Journal of Neurophysiology 38: 356-368)—that amounts to considering the extracellular medium as purely resistive— leads to a more general, time dependent, version of equation 8:
\begin{align}
\label{eq:stat7}\tag{9}
\Phi_e(X,Y,Z,t) = \frac{a^2 \sigma_i}{4 \sigma_e} \int_{x_{min}}^{x_{max}} \frac{1}{\sqrt{(x-X)^2+Y^2+Z^2}}
\frac{\partial^2 V_m(x,t)}{\partial x^2} dx ,.
\end{align}
Notice that the derivation of equations 8 and 9 does not assume anything about the origin of the membrane potential non-uniformity.
If the membrane potential deviation with respect to rest, (\Delta{}V_m), and its derivatives are null at the boundaries of the integration domain, then two rounds of integration by part give (with (X=0) and (h = \sqrt{Y^2+Z^2})):
\begin{align} \label{eq:statPart}\tag{10} \Phi_e(h) = \frac{a^2 \sigma_i}{4 \sigma_e} \int_{x_{min}}^{x_{max}} \left(\frac{3 u^2}{(u^2+h^2)^{5/2}} - \frac{1}{(u^2+h^2)^{3/2}}\right) \Delta{}V_m(u) du , . \end{align}
At that stage, in order to go further, we need an explicit expression or value for (\Delta{}V_m). We are going to solve numerically the H & H equations for that.
We follow here the exposition of Tuckwell (1988) Introduction to theoretical neurobiology. Volume 2. CUP, pp 54-70. We want to solve the following set of equations:
\begin{align} C_m , \frac{\partial V_m}{\partial t} &= \frac{a \sigma_i}{2} \frac{\partial^2 V_m}{\partial x^2} + \overline{g}K n^4 (V_K-V_m) + \overline{g}{Na} m^3 h (V_{Na}-V_m) + g_l (V_l - V_m) + I_A , , \label{eq:HH-PDE}\tag{11}\ \frac{\partial n}{\partial t} &= \alpha_n (1-n) - \beta_n n , , \label{eq:HH-n}\tag{12}\ \frac{\partial m}{\partial t} &= \alpha_m (1-m) - \beta_m m , , \label{eq:HH-m}\tag{13}\ \frac{\partial h}{\partial t} &= \alpha_h (1-h) - \beta_h h , , \label{eq:HH-h}\tag{14} \end{align}
where (C_m) is the membrane capacitance per unit area [F/cm(^2)]; (\overline{g}K), (\overline{g}{Na}) and (g_l) are the potassium, sodium and leak conductances per unit area [S/cm(^2)]; (V_K), (V_{Na}) and (V_l) are the potassium, sodium and leak currents reversal potentials [mV] and (I_a) is "externally" applied current per unit area [mA/cm(^2)].
We will consider a reaction-diffusion system with the form:
\begin{align} \mathbf{u}t = \mathbf{D} , \mathbf{u}{xx} + \mathbf{F}(\mathbf{u}) , , \label{eq:reaction-diffusion}\tag{15} \end{align}
where the (t) subscript stands for the partial derivative with respect to time, the (xx) subscripts stands for the second partial derivative with respect to position, (\mathbf{u} = \left(u_1(x,t),\ldots,u_n(x,t)\right)^T \in \mathbb{R}^n), (\mathbf{D}) is a diagonal (n \times n) matrix of diffusion coefficients (\left(D_1,\ldots,D_n\right)) and (\mathbf{F}(\cdot) = \left(F_1(\cdot),\ldots,F_n(\cdot)\right)^T) is a vector-valued function. The corresponds with the above H & H equations is obtained by setting: (\mathbf{u} = \left(V_m,n,m,h\right)^T); (\left(D_1,D_2,D_3,D_4\right) = \left(\frac{a \sigma_i}{2 C_m},0,0,0\right)), (F_1(\mathbf{u}) = \overline{g}K n^4 (V_K-V_m) + \overline{g}{Na} m^3 h (V_{Na}-V_m) + g_l (V_l - V_m) + I_A), (F_2(\mathbf{u})), (F_3(\mathbf{u})) and (F_4(\mathbf{u})) are given by equations 12, 13 and 14.
Let us consider a simpler problem, the heat equation:
\begin{align} u_t = D , u_{xx} , , \label{eq:heat-equation}\tag{16} \end{align}
where (u(x,t)) is a scalar. A numerical integration procedure is possible by finite differencing. Here, the heat equation (16) is replaced by a finite difference equation whose solution approximates the one of the heat equation. We discretize the (x) axis using (m+1) equally spaced points (with a step (\Delta{}x)) and the (t) axis using (n+1) equally spaced times (with a step (\Delta{}t)). We write the approximate solution as:
\begin{align} U_{i,j} = u(i \Delta{}x,j \Delta{}t) \quad i = 0,\ldots,m ; i = 0,\ldots,n , . \label{eq:discrete-u}\tag{17} \end{align}
The finite difference approximations of the required derivatives are:
\begin{align} u_t(x,t) &\approx \frac{U_{i,j+1}-U_{i,j}}{\Delta{}t} , , \label{eq:u_t}\tag{18} \ u_x(x,t) &\approx \frac{U_{i+1,j}-U_{i,j}}{\Delta{}x} , , \label{eq:u_x}\tag{19} \ u_{xx}(x,t) &\approx \frac{u_x(x,t)-u_x(x-\Delta{}x,t)}{\Delta{}x} , , \nonumber \ &\approx \frac{U_{i+1,j}-2 , U_{i,j} + U_{i-1,j}}{\Delta{}x^2} , . \label{eq:u_xx}\tag{20} \ \end{align}
The numerical integration of the heat equation with the finite difference equation is obtained by establishing a relation between the (U_{i,j+1}) and the (U_{i,j}). One methods approximates the second spatial derivative at (t) by the one at (t+\Delta{}t) giving the scheme:
\begin{align} \frac{U_{i,j+1}-U_{i,j}}{\Delta{}t} = \frac{D}{\Delta{}x^2} \left(U_{i+1,j+1}-2 , U_{i,j+1} + U_{i-1,j+1}\right), . \label{eq:Ames-scheme}\tag{21} \end{align}
Crank and Nicolson use the average of the approximations to the second space derivatives at the (jth) and ((j+1)th) time points to get:
\begin{align} \frac{U_{i,j+1}-U_{i,j}}{\Delta{}t} = \frac{D}{2 \Delta{}x^2} \left(U_{i+1,j+1}-2 , U_{i,j+1} + U_{i-1,j+1} + U_{i+1,j}-2 , U_{i,j} + U_{i-1,j}\right), . \label{eq:Crank-Nicolson}\tag{22} \end{align}
More generally a weight factor (\lambda) can be used with weight (\lambda) for the ((j+1)th) time points and weight ((1-\lambda)) for the (jth) with (0 \le \lambda \le 1). Then with:
\begin{align} r \doteq \frac{D \Delta{}t}{\Delta{}x^2} , , \label{eq:step-ratio}\tag{23} \end{align}
we have:
\begin{align} -r \lambda U_{i-1,j+1} + (1+2 r \lambda) U_{i,j+1} -r \lambda U_{i+1,j+1} = r (1-\lambda) U_{i-1,j} + \left(1-2 r (1-\lambda)\right) U_{i,j} r (1-\lambda) U_{i+1,j}, , \label{eq:general-Crank-Nicolson}\tag{23} \end{align}
where all the unknown terms in (j+1) are on the left side. Since (i = 0,1,\ldots,m) there are (m+1) equations with (m+1) unknown. This integration scheme is called implicit because a linear system must be solved to obtain the values of (u(x,t)) at the next time step. The system defined by equation 23 is tridiagonal and can be solved without matrix inversion. In Python, the scipy.linalg sub-module provides the solve_banded function to work efficiently with linear systems exhibiting a banded structure.