## Numerical Solving of Diffusion Equation with Coordinate-dependent Chemical Potential

Diffusion is the net movement of a substance from a region of higher concentration to a region of lower concentration and is driven by a concentration gradient. The diffusion process is described by the diffusion equation, which provides the coordinate and time-dependent concentration. The diffusion equation in three dimensions reads:

\(\frac{\partial u\left(x,y,z,t\right)}{\partial t}=D(\frac{{\partial }^{2}u\left(x,y,z,t\right)}{{\partial x}^{2}}+\frac{{\partial }^{2}u\left(x,y,z,t\right)}{{\partial y}^{2}}+\frac{{\partial }^{2}u\left(x,y,z,t\right)}{{\partial x}^{2}}\) )

where *u(x,y,z,t)* is the coordinate and time-dependent concentration, and *D* stands for diffusion coefficient. The equation above assumes that the diffusion coefficient is direction and concentration-independent. In complex environments, for example, a nerve cell, additional complexity is added by the coordinate-dependent chemical potential that substantially alters time-dependent concentration profiles for local anaesthetics. Local anaesthetics are lipophilic substances that prefer lipophilic compartments (axolemma, Schwann cell) over aqueous compartments (axoplasm, extracellular fluid). The equilibrium constant is also pH dependent. In addition, because of the cylindrical symmetry of the myelinated axon, it is worth proceeding with cylindrical coordinates. For considering the coordinate-dependent chemical potential, we followed Livshits approach that, in one dimension, reads\(\frac{\partial u\left(x,t\right)}{\partial t}=D\left[\frac{{\partial }^{2}u\left(x,t\right)}{{\partial x}^{2}}+\frac{1}{{k}_{B}T}\left(\frac{\partial u\left(x,t\right)}{\partial t}\bullet \frac{\partial \left(x\right)}{\partial x}+u\left(x,t\right)\bullet \frac{{\partial }^{2}\left(x\right)}{{\partial x}^{2}}\right)\right]\)

where \(u\left(x,t\right)\) is coordinate and time-dependent concentration and \(\left(x\right)\) is coordinate-dependent chemical potential, \({k}_{B}\) is Boltzmann constant, and T is the absolute temperature. The approach is a special case of the Smoluchowski equation for isotropic medium in one dimension.

In cylindrical coordinates, the same equation reads

$$\frac{\partial u\left(R,t\right)}{\partial t}=D\left[\frac{{\partial }^{2}u\left(R,t\right)}{{\partial R}^{2}}+\frac{1}{R}\frac{\partial u\left(R,t\right)}{\partial R}+\frac{1}{{k}_{B}T}\frac{\partial u\left(R,t\right)}{\partial R}\bullet \frac{\partial \left(R\right)}{\partial R}+\frac{1}{{k}_{B}T}u\left(R,t\right)\bullet \frac{{\partial }^{2}\left(R\right)}{{\partial R}^{2}}\right]$$

where R is the radius from the axon centre. Introducing cylindrical coordinates is necessary because the outer layer, consisting of axon membrane and about forty wraps of the Schwann cell, is no longer very small relative to the radius.

Since the diffusion equation with coordinate-dependent chemical potential described above does not allow for an analytical solution, we proceeded with numerical solving on a grid. Algorithms for numerical solving of the diffusion equation with coordinate-dependent chemical potential appeared only recently (19).

The following scheme was applied to integrate concentration as a function of coordinate and time. Fast Fourier transform (FFT) was applied to calculate the first and second derivatives of \(u\left(R,t\right)\) with respect to R at each time step. The first and second derivatives of chemical potential with respect to the coordinate were time-independent and were also calculated with FFT. The applied numerical scheme belongs to the category of spectral methods for solving partial differential equations and is more stable than Crank-Nicolson or forward time-centered space schemes. Spectral methods proved to be very efficient for solving the time-dependent Schrödinger Eq. (20–22). Recently, alternative numerical schemes have been developed for solving diffusion equations with coordinate-dependent chemical potential. Koslowski and coworkers applied a closely related scheme to numerically solve the diffusion equation for cocaine diffusion across the membrane with a Monte Carlo-type approach (23). Ghysels et al. studied oxygen transport across the membranes by solving the Smoluchowski equation, where the potential of mean force was calculated by molecular dynamics simulation in conjunction with biased sampling (24). Please note that in our approach, we are formally changing protonation states of local anaesthetics, which would make the potential of mean force construction from molecular dynamics simulation extremely demanding.

The following propagation scheme was applied

$${u}_{i}^{l+1}={u}_{i}^{l}+D\varDelta t\left[\frac{{\partial }^{2}{u}_{i}}{{\partial R}^{2}}+\frac{1}{R}\frac{{\partial u}_{i}}{\partial R}+\frac{1}{{k}_{B}T}\bullet \frac{{\partial u}_{i}}{\partial R}\bullet \frac{\partial {}_{i}\left(R\right)}{\partial R}+\frac{1}{{k}_{B}T}{u}_{i}\frac{{\partial }^{2}{}_{i}\left(R\right)}{{\partial R}^{2}}\right]$$

where index l = 1…M runs over time and index i = 2…N-1 runs over the coordinate values, and index *i* refers to the values of *u* and \(\)and the corresponding derivatives at the coordinate values corresponding to the index of *i*.

On the left-hand side, we applied the von Neumann boundary condition \(\frac{\partial u\left(x,t\right)}{\partial x}=0\)implying that local anaesthetic flow is zero at the centre of the neuron. Numerically this was implemented as

$${u}_{N}^{l+1}=D\left[{u}_{N}^{l}+\frac{\varDelta t}{{\varDelta x}^{2}}\left({u}_{N-1}^{l}-{u}_{N}^{l}+{u}_{i+1 }^{l}+{h}_{1}\varDelta x\right)\right]$$

where \({h}_{1}=0\).

On the right-hand side of the interval, the Dirichlet boundary condition was implemented, stating that the system is in contact with an infinite reservoir of local anaesthetic solution with the time-independent concentration \({g}_{0}\). Numerically, the Dirichlet boundary condition was implemented with the following equation

$${u}_{1}^{l+1}={g}_{0}$$

We applied a one-dimensional grid with a spacing of 1 Å, which provides adequate precision for a neural cell model. Two numerical experiments were performed. For studying the time to onset, we considered only the membrane at the node of Ranvier (axolemma), which is crucial for signal transduction. For studying the duration of action, we considered the full complexity of the axon membrane and 40 wraps of the Schwann cell.

The applied potential of mean force is shown in Fig. 3

Experimental values for diffusion coefficients of lidocaine and bupivacaine in water of 7.49.1010 Å2 s− 1 and 6.71.1010 Å2 s− 1, respectively, were used (18). The applied time-step for integration was 1 ps which still allows for numerical stability. We integrated the modified diffusion equation for 9 ns, which required less than a day of CPU time. Please note that we considered only the rate of filling of the lipophilic parts from aqueous compartments *k**on*; the corresponding voiding is a much slower process, and its rate constant *k**off* is connected with *k**on* by the principle of detailed balance, where *k**on**/ k**off* *= K**ow*. An equilibrium constant, *K**ow*, was approximated by an experimental octanolwater partition coefficient.

## Construction of a System of Ordinary Differential Equations to Describe Local Anaesthetics’ Micropharmacokinetics

Since the application of a fine grid was required for a realistic description of the myelinated axon, a very small integration step was required. The available computer power did not allow for a long enough simulation time to be relevant for studying the time to onset and duration of action. In order to simulate the time span of hours, we approximated the axon as consisting of four compartments; the axoplasm, node of Ranvier membrane (axolemma), Schwann cell, and extracellular fluid. The local anaesthetic concentrations in the distinct compartments are time-dependent but uniform. The model was described by a system of ordinary differential equations. A scheme of compartments used in the construction of ordinary differential equations, along with the corresponding volumes and labels of permeation rates, is shown in Fig. 4.

We paid particular attention to the microanatomy of peripheral human nerves that are typical targets for local anaesthetics. An example of a micrograph of a transversal section through human *N. ischiadicus* is shown in Fig. 5.

Substantial variability exists between nerves even of the same person; therefore, our model has been based on average literature values. Our model is established on a human nerve with several fascicles. Each fascicle is a collection of axons, Schwann cells surrounding the axons, connective tissue called endoneurium, and capillaries. One such fascicle has an average diameter on the order of magnitude of 350 µm. The fascicle size does not directly influence our model because we model a single axon. Based on the literature data, a typical fascicle has a myelinated axon density of 25969 mm− 2 (25) and a capillary density of 67.9 mm− 2 (26). Furthermore, about 60% of the fascicular cross-section area are axons with Schwann cells; the rest is extracellular fluid (ECF) (25, 27). Because each fascicle encapsulates neurons of various sizes, only small pain-transducing Aδ neurons will be considered, which have an axon with a diameter of 0.5 µm and a myelin sheath with an external diameter of 2 µm. The Aδ neurons are responsible for propagating sharp localized pain in mammals. Since this is the specific pain modality that is usually tested in clinical practice and clinical research on the effectiveness of local anaesthetics, we believed it would be the best model for comparison to the results in the literature. Aδ neurons consist of one continuous axon and many Schwann cells wrapped around the axon, with a small gap between the Schwan cells, called the Node of Ranvier (Fig. 2). We could then consider that the nerve consists of a number of repeatable segments, where each segment consists of a 1501.5 µm segment of axoplasm, 1500 µm long Schwann cell, 1.5µm long Node of Ranvier, and the extracellular fluid, as shown in Fig. 4. Considering this thin segment, we can calculate the cross-sectional surface area of each compartment in the segment followed by the volume of each part of the segment (Table 1).

Table 1

Cross-sectional surface area and volume of parts of the neuron segment consisting of 1501.5 µm segment of axoplasm, one node of Ranvier, one Schwann cell, and the corresponding extracellular fluid.

Compartment | Surface area [µm2] | Volume [µL] |

Axoplasm | 1.963 x 10− 1 | 2.948 x 10− 7 |

Axolemma | 7.933 x 10− 3 | 1.190 x 10− 11 |

Schwann cell | 2.790 | 4.173 x 10− 6 |

Extracellular fluid (ECF) | 1.984 | 2.979 x 10− 6 |

The compartments shown in Fig. 4 are kinetically coupled, and the dynamics are controlled by permeation rate constants. The permeation rate constants have units s− 1 and are scaled for lipophilic compartment volumes (axolemma or Schwann cell). They were obtained by solving the diffusion equation considering experimental geometric parameters of the nerve. The applied estimation of the capillary clearance is based on the published experimental data. The capillary clearance was approximated with the following assumptions: Let the diffusion of local anaesthetic into the capillary be diffusion limited; therefore, the clearance of local anaesthetic by the capillary is equal to the flow rate of the blood in the capillary (28). Typical blood flow velocity in an endoneurial capillary is 0.79 mm s− 1 (29). For a human endoneurial capillary with a lumen of 33 µm2 (26), the flow rate is, therefore, 2.61 x 10− 5 µL s− 1. The ratio between the number of capillaries and the number of neurons is \(\frac{67.9{mm}^{-2}}{25969{mm}^{-2}}=0.0026\). Moreover, normally 25% of capillaries in endoneurium are not perfused at any given time (26). Capillary clearance (*Cl*) per neuron is then calculated as follows:

$$Cl=flowrate\bullet \left(ratio of capillaries to neurons\right)\bullet (1-0.25)=2.61\bullet {10}^{-5}\mu L {s}^{-1}\bullet 0.0026\bullet 0.75=5.09\bullet {10}^{-8}\mu L {s}^{-1}$$

Since the volume of ECF is 2.979 x 10− 6µL, the capillary permeation rate is

$${k}_{c}=\frac{Cl}{V}=\frac{5.09\bullet {10}^{-8}\mu L {s}^{-1}}{2.979\bullet {10}^{-6}\mu L}=0.0171{s}^{-1}$$

We described the time dependence of concentrations in four compartments with a system of four coupled ordinary differential equations, where C1 and V1 correspond to ECF, C2 and V2 correspond to the Schwann cell, C3 and V3 correspond to the axolemma, and C4 and V4 correspond to axoplasm. For a definition of the permeation rate constants, see Fig. 4, while values are collected in Table 2.

$$\frac{d{C}_{1}}{dt}=\left[\left(-{k}_{12}\bullet {C}_{1}\left(t\right)+{k}_{21}\bullet {C}_{2}\left(t\right)\right)\bullet {V}_{2}+\left(-{k}_{13}\bullet {C}_{1}\left(t\right)+{k}_{31}\bullet {C}_{3}\left(t\right)\right)\bullet {V}_{3}\right]\bullet \frac{1}{{V}_{1}}-{k}_{c}\bullet {C}_{1}\left(t\right)$$

$$\frac{d{C}_{2}}{dt}=\left({k}_{12}\bullet {C}_{1}\left(t\right)-{k}_{21}\bullet {C}_{2}\left(t\right)\right)+\left(-{k}_{24}\bullet {C}_{2}\left(t\right)+{k}_{42}\bullet {C}_{4}\left(t\right)\right)$$

$$\frac{d{C}_{3}}{dt}=\left({k}_{13}\bullet {C}_{1}\left(t\right)-{k}_{31}\bullet {C}_{3}\left(t\right)\right)+\left(-{k}_{34}\bullet {C}_{3}\left(t\right)+{k}_{43}\bullet {C}_{4}\left(t\right)\right)$$

$$\frac{d{C}_{4}}{dt}=\left[\left({k}_{24}\bullet {C}_{2}\left(t\right)-{k}_{42}\bullet {C}_{4}\left(t\right)\right)\bullet {V}_{2}+\left(-{k}_{43}\bullet {C}_{4}\left(t\right)+{k}_{34}\bullet {C}_{3}\left(t\right)\right)\bullet {V}_{3}\right]\bullet \frac{1}{{V}_{4}}$$

The system does not allow for an analytical solution; therefore, it was solved numerically using the Runge-Kutta 4th order method, which was done using Fortran 90. The calculations were performed on a desktop computer in double precision. We performed two simulations, one for the absorption phase and one for the voiding phase. The latter mainly depends on the capillary clearance *k**c* in the first differential equation. The absorption phase lasts until the equilibrium is reached. The time step for numerical integration was 0.1 µs when studying the absorption phase and 0.1 ms when studying the voiding phase. The integration covered the real-time span of hours when studying the voiding phase. For the simulation of the voiding phase, we disregarded the axolemma compartment, allowing for a larger time step. This was done by setting the permeation rates for this compartment to be zero. The described simplification did not affect the time course of concentration in any other compartment during the voiding phase because the axolemma compartment is 10000 times smaller than other compartments; therefore, its storage capacity can be neglected, which we proved numerically. Nevertheless, the short time step was necessary because of the first few seconds of integration. Since the calculations took only a couple of hours of CPU time, we did not improve the applied algorithm with adaptive time step integrators.