PyPNS: Multiscale Simulation of a Peripheral Nerve in Python

Abstract

Bioelectronic Medicines that modulate the activity patterns on peripheral nerves have promise as a new way of treating diverse medical conditions from epilepsy to rheumatism. Progress in the field builds upon time consuming and expensive experiments in living organisms. To reduce experimentation load and allow for a faster, more detailed analysis of peripheral nerve stimulation and recording, computational models incorporating experimental insights will be of great help. We present a peripheral nerve simulator that combines biophysical axon models and numerically solved and idealised extracellular space models in one environment. We modelled the extracellular space as a three-dimensional resistive continuum governed by the electro-quasistatic approximation of the Maxwell equations. Potential distributions were precomputed in finite element models for different media (homogeneous, nerve in saline, nerve in cuff) and imported into our simulator. Axons, on the other hand, were modelled more abstractly as one-dimensional chains of compartments. Unmyelinated fibres were based on the Hodgkin-Huxley model; for myelinated fibres, we adapted the model proposed by McIntyre et al. in 2002 to smaller diameters. To obtain realistic axon shapes, an iterative algorithm positioned fibres along the nerve with a variable tortuosity fit to imaged trajectories. We validated our model with data from the stimulated rat vagus nerve. Simulation results predicted that tortuosity alters recorded signal shapes and increases stimulation thresholds. The model we developed can easily be adapted to different nerves, and may be of use for Bioelectronic Medicine research in the future.

Keywords: Bioelectronic medicines; Biophysics; Finite element model; Peripheral nerve; Simulation.

Fig. 1

The validation data were obtained through stimulation of a rat vagus nerve. A pseudotripolar electrode excited axons at the cervical vagus nerve, signals were picked up at the subdiaphragmatic vagus nerve with a bipolar electrode

Fig. 2

The Axon-class is the central object of PyPNS’s internal information flow. Together with its associated ExcitationMechanism s it defines the NEURON simulation. Extracellular-objects allow the calculation of extracellular potentials given current i(t), source position s and receiver position r. They are used by both StimField for extracellular stimulation and by RecordingMechanism for recording. All classes are managed in the Bundle-class and supported by helper modules spikeTrainGeneration, signalGeneration and createGeometry.

Fig. 3

Linear and quadratic fits were used to extrapolate the parameters of myelinated axons to smaller diameters. a Diameters of all segments – nodes, MYSA (myelin attachment segments), and paranodal elements FLUT (paranode main segment) and STIN (internode segment, see McIntyre et al. (2002) for more information on the model)–were fit quadratically to prevent negative values. Node distance (b) and number of myelin sheaths (c) were extrapolated linearly

Fig. 4

Axon segments can be interpreted as current point sources. The extracellular potential ϕ(r,t) at position r caused by a current i(s,t) at position s is determined by current time course scaled with a static potential depending on the extracellular space and the spatial relation between source and receiver position

Fig. 5

The impact of the longitudinal profile ϕ_SFAP(z) on SFAPs can be understood by studying the potential caused by a perfectly straight axon recorded at z₀ = 0 for $t=t^{\prime }$. Axon segments of length Δz exhibit the exact same current time course except for a delay Δt = z_i/CV (a). The potential ϕ_SFAP at $t=t^{\prime }$ is then obtained as the sum over membrane currents $i\left(z_i\right|t=t^{\prime })$ shown in (b), multiplied by the static potential ϕ_static(z_i,I_ref)/I_ref (c)

Fig. 6

A circularly symmetric geometry makes it possible to import precomputed potential fields. The nerve is modelled as axons (white matter) surrounded by the epineurium. The positions of exemplary current point sources, each generating one potential field, are shown. For radially inhomogeneous media, a line of sources does characterise all unique fields. For longitudinal inhomogeneities (a), potential fields for a two-dimensional array of point current sources need to be precomputed (b)

Fig. 7

Unmyelinated axons (a) produce a smoother membrane time course than myelinated (b) ones. Both axons had a diameter of 3 μ m. c Unmyelinated axons produce a higher current output per distance. The integrated absolute current during a single action potential over axon length is shown

Fig. 8

Compared to the homogeneous and radially inhomogeneous extracellular media the cuff insulation caused a much softer and strikingly linear characteristic. Electrode radius was 235 μ m. The profiles are shown for two radial axon displacements in solid and dashed lines respectively

Fig. 9

Unmyelinated and myelinated SFAPs showed different sensitivities towards the extracellular space. In the upper plots (a, b), diameters were set to 3 μ m. a The main peak of unmyelinated fibres mostly varied in amplitude over media, not in shape. In cuff insulated nerves, additional side peaks emerged. b Myelinated fibres produced much higher and longer lasting SFAPs in the cuff insulated medium. Both axons had diameter 3 μ m, were placed centrally within the nerve and recorded by a circular monopolar electrode with radius 235 μ m. Conductivity of the homogeneous medium was set to 1 S m^− 1. Lower row shows zoomed-in plots. c The amplitude boost achieved by cuff insulation was stronger for myelinated than for unmyelinated axons over the whole diameter range. For the other two media, unmyelinated SFAPs produced stronger SFAP amplitudes at diameters above 0.5 and 1 μ m respectively

Fig. 10

Unmyelinated and myelinated axons have different ideal cuff lengths. For unmyelinated fibres (a), very short ranges around 1 mm produce the maximal amplitude. For myelinated ones (b), the amplitude only rises with length. Contour lines show the peak-to-peak amplitude in μ V

Fig. 11

a The simulated compound action potential in the cuff medium approaches the experimental recording well in the relevant signal segments. As expected, homogeneous and radially inhomogeneous media lead to much weaker signal amplitudes. For the experimental recording, the grey underlying area indicates the standard deviation over the 10 stimulation repetitions. See Table 2 for axon properties. Distance between stimulation site and bipolar electrode (3 mm pole distance, 235 μ m radius) was 8 cm. All axons were activated by intracellular stimulation. The timing of unmyelinated SFAPs was adapted to regular conduction velocity values assumed in mammalian peripheral nerves (CV $=1.4\cdot \sqrt{d}$, CV in ms^− 1, d in μ m). b The signal from myelinated fibres, which arrive first, appears similar to the experiment. c The unmyelinated signal segment matches the amplitude and duration of the experimental recording as well. The signal-to-noise ratio of the recordings is much worse for unmyelinated fibres, however, as the amplitude of their SFAPs is low

Fig. 12

In the frequency domain, simulation and experiment did not match equally well for both fibre types. a For unmyelinated axons, the simulation did not perfectly approach the experimental spectrum in any medium with best results for the cuff. b The simulated frequency characteristic of myelinated axons in the cuff insulated medium was close to reality

Fig. 13

Fluorescence microscopy images of the mouse sciatic and vagus nerve both show slight tortuosity in their axon trajectories. a The thick myelinated fibres in the sciatic nerve appear very parallel. b The thinner axons in the vagus take a more curvy trajectory. Several manually traced fibres used to fit the model are highlighted in orange

Fig. 14

The axon placing algorithm result (b, c) was fit to tortuosity of microscopy imaged fibres (a). a Direction change distributions (c-distributions) for vagus and sciatic nerve. b c-distributions in the simulation for both normally and uniformly distributed amplitude of the random component ||w|| in Eq. 3 for α s of 0.2, 0.6, and 1.0. c Example axon trajectories in space for uniform (upper) and Gaussian (lower) ||w||-distributions at α-values of 0.2, 0.6, and 1.0

Fig. 15

Unmyelinated axons were more sensitive to tortuousity in their SFAP shapes than myelinated ones. Tortuosity parameter α was set to 0.2, 0.6, and 1.0 for the signals shown in the first two rows. Grey lines correspond to SFAPs of different trials (axon geometries) at the same parameters. a Unmyelinated axons produced SFAPs differing both in timing and shape for the non-insulated nerve (radially inhomogeneous medium), even for small α-values of 0.2. In the cuff-insulated nerve (middle row) their signals became noisy at low α-values and the main SFAP peak almost disappeared for α = 1.0. b Myelinated axons mostly differed in timing in the radially inhomogeneous extracellular space, and not as much in shape. In the cuff, noisiness only arose at high tortuosity values. In the lower plots, the mean maximum pairwise cross-correlation gives a quantitative confirmation of the higher susceptibility of unmyelinated axons to change their SFAP shape in the presence of tortuosity. Note the different ordinate scales

Fig. 16

Unmyelinated axons have higher stimulation thresholds and are activated less reliably than myelinated ones. Both bundles consisted of 15 axons with diameter 3 μ m and were stimulated with a bipolar electrode of radius 235 μ m and pole distance 1 mm using a biphasic pulse of frequency 1 kHz, duration 1 ms and duty cycle 0.5. The extracellular medium was a nerve of diameter 240 μ m bathed in oil. a Unmyelinated axons started to be activated at 1 mA and showed a peak in activation ratio at about 3 mA. b Myelinated fibres had a sharp activation threshold at a much lower current of about 0.15 mA and stayed activated for higher currents. Only when incrementing the stimulation current in very small steps of about 10 nA c a slight tortuosity-induced increase in stimulation threshold became visible for them as well

Fig. 17

If bundle b_k and current axon segment a_i have a fixed relation, e.g. parallel, the expected distribution of segment direction differences c = ||a_i −a_i+ 1|| can be easily obtained from the distribution of ||w|| ($\mathcal{P}$) by their geometrical relation

Fig. 18

An analytic transmission function implements the relation between current source position and potential for recording in a cuff electrode. In the shown case, a nerve of diameter 480 μ m in a cuff of 2 cm length was simulated in the FEM model. Functions are displayed for three different angles between the perpendiculars of source and electrode position onto the bundle guide respectively