Tutorial: using NEURON for neuromechanical simulations

The dynamical properties of the brain and the dynamics of the body strongly influence one another. Their interaction generates complex adaptive behavior. While a wide variety of simulation tools exist for neural dynamics or biomechanics separately, there are few options for integrated brain-body modeling. Here, we provide a tutorial to demonstrate how the widely-used NEURON simulation platform can support integrated neuromechanical modeling. As a first step toward incorporating biomechanics into a NEURON simulation, we provide a framework for integrating inputs from a “periphery” and outputs to that periphery. In other words, “body” dynamics are driven in part by “brain” variables, such as voltages or firing rates, and “brain” dynamics are influenced by “body” variables through sensory feedback. To couple the “brain” and “body” components, we use NEURON's pointer construct to share information between “brain” and “body” modules. This approach allows separate specification of brain and body dynamics and code reuse. Though simple in concept, the use of pointers can be challenging due to a complicated syntax and several different programming options. In this paper, we present five different computational models, with increasing levels of complexity, to demonstrate the concepts of code modularity using pointers and the integration of neural and biomechanical modeling within NEURON. The models include: (1) a neuromuscular model of calcium dynamics and muscle force, (2) a neuromechanical, closed-loop model of a half-center oscillator coupled to a rudimentary motor system, (3) a closed-loop model of neural control for respiration, (4) a pedagogical model of a non-smooth “brain/body” system, and (5) a closed-loop model of feeding behavior in the sea hare Aplysia californica that incorporates biologically-motivated non-smooth dynamics. This tutorial illustrates how NEURON can be integrated with a broad range of neuromechanical models. Code available at https://github.com/fietkiewicz/PointerBuilder.

conceptually distinct biomechanical elements to interact in a natural way. 23 In addition to the conceptual reorganization of the code, another difference is that the Kim model used 24 a multicompartmental, spatially reconstructed motoneuron with customized ion channel mechanisms, 25 whereas the model presented here uses the hh mechanism which is part of the standard NEURON library.  The equations for calcium dynamics and muscle activation are as in Kim (Kim, 2017). We reproduce the 31 equations below, and provide descriptions of the variables in Table 1. Table 2 provides descriptions of the 32 parameter values. 33Ċ a SR = −K 1 · CS 0 · Ca SR + (K 1 · Ca SR + K 2 ) · Ca SRCS − R + U (1) Ca SP = −(K 3 · B 0 + K 5 · T 0 ) · Ca SP + (K 3 · Ca SP + K 4 ) · Ca SPB + (3) Ca SPB = K 3 · B 0 · Ca SP − (K 3 · Ca SP + K 4 ) · Ca SPB (5) Ċ a SPT = K 5 · T 0 · Ca SP − (K 5 · Ca SP + K 6 ) · Ca SPT (6) U = U max · (Ca SP ) 2 · K 2 1 + Ca SP · K + (Ca SP ) 2 · K 2 2 (8) ϕ (X m ) = ϕ 1 · X m + ϕ 2 , for X m < −8 ϕ 3 · X m + ϕ 4 , for X m ≥ −8 (10) A ∞ = 0.5 1 + tanh  It is common to begin with an existing mod file that has a system of equations but does not use pointers.

72
Here we present a particularly simple example of converting a single mod file into multiple files through 73 the use of pointers. The model consists of the two-dimensional, nonlinear system of the Lotka-Volterra 74 equations, also known as the predator-prey equations, given below: where a is interpreted as a prey species, and b is interpreted as a predator species. For concreteness, we 76 will use parameter values α = 1.1 (prey intrinsic birth rate), β = 0.4 (predation effect on prey), γ = 0.1 77 (predation benefit to predator), δ = 0.4 (predator intrinsic death rate).

78
Using a single file, the model above might be implemented as shown below, where the mechanism is 79 named model.  Next, we place each differential equation in a separate file using a POINTER statement in the NEURON 102 block and placing the name of the pointer variable in the ASSIGNED block. Figure 1 shows the 103 implementation using separate mechanisms named prey and predator.  Writing a setpointer instruction has two particular complications. First, it requires the names of multiple 122 variables, mechanisms, and sections, for a total of six separate pieces of information. Second, the syntax 123 for combining these pieces takes multiple forms, depending on the type of variable (state or parameter) and 124 programming language (hoc or Python). We present two GUI applications, for hoc and Python respectively, 125 that can help a modeler organize the required components and automatically generate the correct syntax.

126
Source code for the tools is provided, and interested users can customize them if desired.

127
The first tool is intended for users who wish to work with the hoc language. It is written using the native 128 NEURON library for GUI development and does not require Python. The interface is shown in Figure 3,  Buttons in the bottom row of the main GUI are used to either save or display the resulting hoc instruction 137 (see Figure 4). The option to "Make .hoc file" saves the hoc instruction to a file which can be automatically loaded by other hoc programs. The option to "View .hoc command" displays the hoc instruction in the 139 window and in the NEURON terminal window. Note that NEURON allows direct text copying from the 140 terminal window but not from the GUI.

141
The second tool is intended for users who wish to work with the Python language. It is written using the 142 native GUI tool for Python, called Tkinter, and is shown in Figure 5. The layout and action buttons are 143 similar to those in the hoc tool (see Figure 3), except that a Python command is generated instead. When 144 viewing the Python command, as seen in Figure 6, the user can copy text directly from the window. The

145
Python command can also be saved to a Python file that can be imported by other Python programs. The example shown is based on the NMODL programs in Figure 1 for the Lotka-Volterra model. The interface is comprised of three groups of controls: Pointer Settings, State/Parameter Settings, and Actions. The "Pointer Settings" indicate that the predator mechanism contains the pointer variable aPointer, and is inserted in the section model. The "State/Parameter Settings" indicate that the prey mechanism contains the state variable a, and is inserted in the section model. Available "Actions" include viewing the hoc command, saving the command as a hoc file, saving and loading Pointer Builder settings, and reading help documentation.   Figure 1 for the Lotka-Volterra model. The interface is comprised of three groups of controls: Pointer Settings, State/Parameter Settings, and Actions. The "Pointer Settings" indicate that the predator mechanism contains the pointer variable aPointer, and is inserted in the section model. The "State/Parameter Settings" indicate that the prey mechanism contains the state variable a, and is inserted in the section model. Available "Actions" include viewing the Python command, saving the command as a Python file, saving and loading Pointer Builder settings, and reading help documentation. Figure 6. Python Pointer Builder output. A complete command for setting the pointer is displayed after selecting the "Python command" option. The text can be copied directly from the window.

147
Section 2.2 presented a neuromechanical, closed-loop model of a half-center oscillator coupled to a rudimentary motor system. For the custom brain and body mechanisms, partial details were provided in section 2.2, and the complete source files are available using the link given in Methods, section 4. For the full model specification, see (Yu and Thomas, 2021). We reproduce the basic equations here. The model comprises two Morris-Lecar type models each with membrane potential V i and potassium gating variable where C is the capacitance, I ext is an applied current, g L is the passive leak conductance, g Ca is the maximal calcium conductance, M ∞ is the voltage-dependent fraction of open calcium channels, E X , X ∈ {Ca, K, L} is the reversal potential for each current, respectively, and g K is the maximal potassium conductance. The coupling function S CPG ∞ (V j ) is a saturating sigmoidal function that specifies the impact of the voltage of cell j on the membrane dynamics of cell i. The strength of the inhibitory synapse forming the HCO is controlled by g CPG syn , and the sensory feedback from the "body" to the "brain" is given by g FB syn S FB ∞ (L j )(V i − E FB syn ). In addition, The force generated by muscle i is Each mucle has a length-tension curve given by where ℓ is a fixed, nominal length value, and −50 < x < 50 is the position of the "limb" driven by the 149 opposing muscle forces. Muscle activation a i is described by the variable where Finally, the position of the "limb" is determined by the balance of forces: The CPG comprises a membrane potential V , and dynamical gating variables n (a delayed rectifier potassium (I K ) channel activation) and h (persistent sodium (I NaP ) channel inactivation). The model sets the "instantaneous" gating variables p ∞ (I NaP activation) and m ∞ (fast sodium (I Na ) activation) to be equal to their voltage-dependent asymptotic values. The I Na inactivation gate is set equal to (1 − n). The model also includes a leak current (I L ) and a tonic excitatory (I tonic ) current. All together, we have

Closed-loop respiratory control model
with parameters C = 21 pF, g K = 11.2 nS, g NaP = 2.8 nS, g Na = 28 nS, g L = 2.8 nS, E K = −85 mV, Motor pool activity: The BRS cell's membrane potential (V ) provides the output of the CPG, driving the respiratory muscles through synaptic activation of a motor unit (α): .
Here, we set r a = r d = 0.001 mM −1 ms −1 giving the rise and decay rate of the synaptic conductance.

163
Lung volume: The output of the motor unit determines the rise and fall of lung volume (vol L ): Here vol 0 = 2 L is the volume of the unloaded lung, and parameters E 1 = 0.4 L and E 2 = 0.0025 ms −1 165 were chosen so that the lung expansion would remain in a physiologically reasonable range (West, 2008).

167
Lung oxygen: At standard atmospheric pressure (760 mmHg), external air with 21% oxygen content registers a partial pressure of oxygen of P ext O 2 = 149.7 mmHg. As the lungs expand d dt [vol L ] > 0 , they draw in external air. The model makes the simplifying assumption that this fresh air mixes instantaneously with the air already present in the lungs. Therefore, the partial pressure of oxygen in the lung alveoli (P A O 2 ) increases at a rate given by the pressure difference between external and internal air, and by the lung volume. During exhalation, on the other hand, d dt [vol L ] ≤ 0 , no external air enters the lungs, so the mixing of air stops. During both contraction and expansion of the lung, oxygen moves between the lungs and the blood. The flux of oxygen from the lungs to the blood occurs at a rate determined by the time constant τ LB = 500 ms, and by the difference in partial pressure of O 2 between the lungs (P A O 2 ) and the arterial blood (P a O 2 ). The rate of change in P a O 2 is given by: where the notation [x] + indicates max(x, 0).

169
Blood oxygen: To represent the change in P a O 2 , we have (The fluxes of oxygen from the blood to the tissues (J LB ) and from the lungs to the blood (J LB ) have units of moles of O 2 per millisecond.) The denominator converts changes in the number of moles of O 2 in the blood to changes in P a O 2 . To calculate the flux J LB , the model assumes the ideal gas law P V = nRT , where n is the number of moles of O 2 , R = 62.364 L mmHg K −1 mol −1 is the universal gas constant, and T = 310 K is temperature. The resulting flux depends on the difference in oxygen partial pressure between the lungs and the blood: Note that the term J BT accounts for both dissolved and hemoglobin-bound oxygen in the blood: Following Henry's law, the model takes the concentration of dissolved oxygen in the blood to be directly proportional to P a O 2 . The blood solubility coefficient, β O 2 = 0.03 ml O 2 × L blood −1 mmHg −1 for blood at 37 degrees C, is the constant of proportionality. Most of the blood's stored oxygen is bound to hemoglobin (Hb). Cooperative binding of oxygen to the four binding sites in each hemoglobin molecule leads to a sigmoidal hemoglobin saturation curve SaO 2 : The phenomenological parameters K = 26 mmHg and c = 2.5 are taken from Keener and Sneyd (2009). connecting P a O 2 with the conductance representing external drive to the CPG (g tonic ): 178 Here, ϕ = 0.3 nS, θ g = 85 mmHg, and σ g = 30 mmHg. This conductance closes the control loop, since The dynamics of the nominal Aplysia feeding model, also known as the SLG model (Shaw et al., 2015; The variables a i represent the firing rates of three nominal populations of motor neurons, active during 189 the "protraction-open" phase (a 0 ), the "protraction-closed" phase (a 1 ), and the "retraction" phase (a 2 ). We (r = 1). This Boolean switching term appears in the differential equation for the grasper position, dx r /dt.

197
The protraction (dx r /dt > 0) or retraction (dx r /dt < 0) of the grasper is determined by the balance of 198 forces produced by the protractor muscle activation (u 0 ) and retractor muscle activation (u 1 ). The net force 199 exerted by the muscles is given by the sum of the two muscle forces 200 F musc (u 0 , u 1 , x r ) = F musc,pro + F musc,ret x(x − 1)(x + 1) is the effective length-tension curve for muscle forces, c i , w i and k i denote the mechanical properties of 202 each muscle.

203
The seaweed is subject to an external force resisting ingestion, F sw . This force only affects the grasper 204 when the grasper is closed (r = 1). This condition is set by a threshold, (r = 1 ⇐⇒ a 1 + a 2 > 0.5).
When the grasper is open (a 1 + a 2 ≤ 0.5), r = 0, and the grasper moves independently of the seaweed.

206
Values for model parameters and initial conditions are given in Table 3 and Table 4   State variable Initial value Description a 0 0.9 activity of I2 motor pool (non-negative) a 1 0.08355 activity of hinge motor pool (non-negative) a 2 0.00003 activity of I3 motor pool (non-negative) u 0 0.748 activity of I2 muscle u 1 0.25 activity of I3 muscle x r 0.65 grasper position (0 is retracted, 1 is protracted) x sw 0 seaweed position (positive is away from the animal)