Incorporation of the MRG model into FEM, units

The basics of how to develop, test, and use models.
Post Reply
Posts: 1
Joined: Thu Apr 13, 2017 10:50 am

Incorporation of the MRG model into FEM, units

Post by ilya747 »


I am trying to incorporate the MRG model (McIntyre, Richardson and Grill, 2002) into COMSOL environment (with FEM) to see the movement of the AP along the fiber represented by a line.

The model does not work properly for me; I think the reason may be in converting the cable equations derived with Kirchhoff' and Ohm' laws into ones operating with "densities" of parameters (R, C, etc).

The equations I've got in absolute values are below. No stimulation current is applied.

1. Axolemma
cm*dv/dt = (1/r_a)*d2v/dx2 - I_ionch - for node, v_my at node = 0;
ca*dv/dt = (1/r_a)*d2(v+v_my)/dx^2 - g_pas*(v - e_pas) - for internode

2. Periaxonal space
xc*d(v_my)/dt = (1/r_a)*d2(v+v_my)/dx2 + (1/r_pa)*d^2(v_my)/dx2 - v_my*xg;

cm, ca, xc - capacitances of node, internodal axolemma and myelin [F/cm];
v, v_my - potential across the axolemma and myelin sheath [mV];
r_a - axoplasmic resistance [kOhm/cm];
g_pas - mysa, flut or stin conductance [mS/cm];
r_pa - resistance of periaxonal space [kOhm/cm].
xg - conductance of myelin [mS/cm]

This is how I convert the parameters to densities provided in the paper:

cm = (S_node/L_node)*Cm[F/cm^2]; S_node = 2*pi*(r_node)*L_node;
ca = (S_ax/L_ax)*Ca[F/cm^2]; S = 2*pi*(r_ax)*L_ax; L_ax = 2*L_mysa + 2*L_flut + 6*L_stin;
xc = (S_my/L_ax)*C_my[F/cm^2]/(2*N); S_my = 2*pi*L_ax*N*(r_ax + r_fib);
Myelin area was found as sum surface area of N cylinders with evenly spaced radii from r_ax to r_fib, N - number of myelin lamellae;
r_a = (L_segm/V_segm)*Ra[kOhm*cm]; V_segm = pi*(r_segm)^2*L_segm;
r_pa = (L_segm/V_segm)*Rpa[kOhm*cm]; V_segm = pi*((r_ax+space)^2-(r_ax)^2)*L_segm;
xg = (S_my/L_ax)*g_my[mS/cm^2]/(2*N); S_my = 2*pi*L_ax*N*(r_ax + r_fib);

If I put all this into the equations applied to 1D FEM model of the fiber, the AP doesn't propagate (just passive fall of the potential is seen), unless I significantly change some parameters, e.g. decrease values of ca and xg. Somewhere must be a mistake.

The parameters in the model I've found in the NEURON db ( ... model=3810) seem to be defined relative to the fiber diameter (with myelin sheath) and not the axon diameter. That is how I understand why there are expressions like "Ra=rhoa*(1/(paraD2/fiberD)^2)/10000".
So, I redefined everything relative to the diameter of the fiber with myelin sheath. The model now stopped converging when the potential reaches threshold (backward differentiation with up to 5th order is used),a mistake must still be there. Size of the mesh does not affect the result as well.

Sorry for so much text in here, I've already spent considerable amount of time on this model and would very much appreciate any help.

Thank you !
Post Reply