Influence of cellular structures of skin on fiber activation thresholds and computation cost

Electrical neuromodulation is widely used to treat and manage neurological disorders. Migraine, a socioeconomic burden, may be treated using this technique. Transcutaneous stimulation of frontal nerves by electrodes placed on the forehead is of interest as it exposes patients to lower levels of risk and side-effects compared with surgical and pharmaceutical solutions and may be readily delivered. The size, shape and placement of the electrodes can be optimised using computational models involving a volume conductor model of anatomical structures and electrodes as well as nerve fibre models. A detailed volume conductor incorporating cell level structures of skin can yield an accurate map of electrical potential distribution due to an electrode setting. However, such a model imposes a very substantial computational cost which may impede the design process. Computation cost can be significantly reduced if the skin microscopic structures are ignored. In this study, we compare the accuracy and computation cost with and without skin microscopic structures on the outcome of a device for transcutaneous frontal nerve stimulation. The performance is presented as the percentage activation of target nerve fibres in response to the level of stimulus current delivered via surface electrodes placed on the forehead. When cell level structures of skin are not incorporated, discretisation time is reduced from 21 h to 0.4 h and the number of finite elements used from 18 M to 1.4 M. Only 1% difference in stimulus current thresholds is observed.


Introduction
Migraine is a highly prevalent and disabling neurological disorder. It is characterized by headache, nausea, vomiting, photophobia and phonophobia [1]. It has a significant impact on the quality of life for those affected [2]. The common complaint of migraine sufferers is generally the symptoms of pain originating in the frontal region of the head [3]. Thus, migraine may be primarily associated with the supraorbital nerve (SON) and supratrochlear nerve (STN) which are branches of the frontal nerve. A solution targeting this specific nerve may be of interest.
Available pharmaceutical and surgical treatments of migraine are generally not completely efficacious and have troublesome side-effects [4]. Trancutaneous neuromodulation, in which current is applied through the skin via surface electrodes to stimulate the central or peripheral neural tissue, has been used for pain management. Among different transcutaneous electrical nerve stimulation (TENS) methods, transcutaneous frontal nerve stimulation (t-FNS) has been applied on a large group of patients who have episodic migraine using a device called Cefaly (Cefaly, CEFALY Technology, Liège, Belgium). By applying stimulus current to sensory fibers in frontal nerve, migraine has been shown to be alleviated [5].
The effectiveness of t-FNS in episodic migraine was tested in a randomised double blind controlled study trial supported by large scale post marketing survey. In nearly half of the patients t-FNS was not efficacious [6,7]. This inconclusive response may be associated with neuro-anatomical variations in patients [8] or electrode arrangements (size and shape). No study has investigated the underlying causes of inefficacies. This is partly due to the physical limitations of studying the neuroanatomy of individuals and different settings of electrode arangements. Computational models [9] may enable researchers to estimate current stimulation thresholds in neuromodulation therapy and investigate the effects of various parameters. Such models are usually implemented in finite element (FE) models (FEM) involving a volume conductor model representing various anatomical structures and the electrodes by their respective conductivities and appropriate boundary conditions. Since the electrode patch, as the interface for delivering the stimulus current, is in contact with the skin, incorporating microscopic details of this layer may be important in simulating the resulting electrical potential fields accurately. Mammalian skin is comprised of epidermis and dermis layers. The outer layer of the epidermis is called stratum corneum (SC) [10]. The thickness of SC is between 10 μm and 50 μm and contains 10-20 layers of dead keratinized cells (keratinocytes) with lipid lamellae filling the intercellular regions and sweat ducts (SDs) [11]. Since the size of these layers are considerably smaller than other tissue layers involved, a large dimensional ratio between these structures will lead to a significant computation cost in FE models. Microscopic structures may be ignored at the cost of a degree of error to provide managable computational requirements.
The effect of cell level modeling of the skin on the electrical field in the volume conductor (VC) and temperature elevation around the electrode have been investigated [12]. In our previous study [13], we compared a realistic human head model based on magnetic resonance imaging (MRI) scans and a geometrically simplified human head model (generated from geometric shapes (e.g., sphere)) with respect to the estimates of stimulus current thresholds. The results show that the simplified model can be used intead of the realistic model with some error (3%, based on stimulus current thresholds). The cellular structures of the skin cannot be evaluated in a simplified model of the entire head due to the computation cost. Therefore, a region of interest (forehead) with all subsequent tissue layers were considered for this study.
In this study, two models were developed for further analysis. In one, the cellular structures of skin were incorporated, referred to as the cellular model, and in the other, the microscopic structures of skin were excluded, here referred to as the simplified model. In the cellular model, the skin was delaminated into multiple layers (as shown in figures 1(b) and (c)). The effect of the lipid lamellae, keratinocytes and uniform random variations of the SDs on the estimated current thresholds and current density of the nerve fibers were studied using hybrid (coupling the FEM results with the Neuron model) computational modeling. The simplified and cell level models were further compared based on electrode size with respect to their current thresholds and computational features to obtain an efficient and reliable trancutaneous electrical stimulation (TES) VC for further studies.
For all the subsequent simulations and operations, a computer with an Intel Core i7-6700 CPU @ 3.4 GHz with 64 GB RAM was used.

Design of the simplified volume conductor
The head tissue layers were built from spherical and TES electrodes were constructed using many alternative geometric shapes in COMSOL (COMSOL, Ltd, Cambridge, UK). The block diagrams of these layers are shown as in figure 1(b). The curvature of the region of interest (forehead) was constructed to follow that of the realistic human head model. To increase model accuracy, the simplified model was constructed based on the average thicknesses of those layers in the MRIbased realistic model (skin: 2.8 mm, fat: 2 mm, muscle: 1.7 mm, skull: 5.5 mm and Cerebrospinal fluid (CSF): 1.5 mm). The white and grey matters were combined and modeled as brain as electrical potential field decays sharply inside them. The trajectory of frontal nerve passes below the frontal bone and exits from the orbital rim and penetrates the corrugator and frontalis muscle. The realistic trajectory of the nerve, derived from anatomical data, was used in both cellular and simple VCs [13]. The SC layer is comparatively thin and was defined using a boundary condition at the outermost boundary of skin to decrease computation cost during simulations in COMSOL.

Design of the cellular volume conductor
The stimulating current flows through the skin via lipid lamellae, keratinocytes and SDs [11]. These layers may have an impact on the potential distributions in the VCs which may affect percentage activation (PA) of the nerve fibers. Therefore, the cell level model was designed to examine the impact of these layers on TES modeling. The block diagram of these layers are shown in figure 1(c).
Stimulation currents generally have the highest intensity underneath electrodes. Thus, only the cellular structures of the skin layer immediately below them were considered. The rest of the tissue layers (e.g., fatty tissue) were of the same dimensions as those in the simplified model. The cellular structures in the skin layer were derived from their typical morphological parameters. However, the values of typical thickness ( δ tt ) of the lipid and keratinocytes were limited by the memory size of the existing PC. The typical and proposed morphological parameters of the cellular structures are detailed in table 1.
The lipid-keratinocytes layers are crossed by SDs as shown in figure 1(c). These SDs start from dermis layer and extend upwards to the surface. The average density of SDs on the forehead of human is 3/mm 2 , their length is approximately 0.3 mm and their diameter is 98±11 μm [15], [16]. Based on the given data, the distributions of the SDs across the human forehead were generated using a uniform random function in MATLAB v.R2015b (MathWorks, Inc., Natic M, USA). SDs with various diameters were generated from smooth geometric shapes (to decrease complexity) and designed accordingly in the VC based on their average density per mm 2 in the region of interest. During modelling of SDs in the cellular model, random spacing was chosen for x and y direction to obtain more realistic results. It is vital to have the same trajectory of the nerve for both models for a fair comparison.
The conductivity of tissue layers were selected based on low frequencies. After generating VCs, the associated conductivity of each tissue layer was obtained to solve underlying equations. The anisotropy of the muscles and SDs were considered in their conductivities. Remaining tissue layers were defined as isotropic. The conductivities of different components is summarized in table 2. Since the skin layer was divided into multiple layers in the cellular model, the conductivities in table 2.a were used for cellular structures. The conductivities of the remaining tissue layers were used as in table 2.b.

Electrode settings
The electrode size was assumed to be smaller than the Cefaly electrode due to the limitations of the computation capabilities of available workstations. Although the same electrode size was used for both models, it is important to investigate the range of electrode sizes that minimally affect the outcome.
It is assumed that the activation of the nerve fibers mostly depends on the extracellular potential variations across the trajectory of the nerve fibers. Thus, if an electrode setting leads to a similar extracellular potentials and PAs versus stimulus current levels variations for both simplified and cellular models, this electrode setting can be used to compare the computation features for both models.
The smallest electrode size was derived from Cefaly electrodes. The parameters of this electrode were proportionally increased until they resulted in  During the electrode design, the nerve trajectory was positioned under the center of the electrode to obtain the maximum possible difference for both models. This procedure was then repeated for each different electrode size. In this study, the largest electrode 'E1' had the smallest error between cellular and simplified models with respect to the electrical potential field and required current thresholds. This size was used to evaluate the current density and computation cost.

Finite element method (FEM) simulation
The domains in the VCs were discretized using free tetrahedral elements to solve underlying differential equations. The region of interest was more finely meshed, while the rest of the region was relatively coarsely meshed to obtain potential distributions in a reasonable time. The total number of elements in each model is compared in figure 3(b). The triangular and edge elements for both simplified and cellular models were 0.085 M to 4 M and 0.005 to ∼1 M, respectively. The average tetrahedral mesh quality and average growth rate for simplified and cellular models were 0.67 to 0.65 and 2.02 to 2.06, respectively. It is noted that the triangular mesh quality was about 0.96 and nearly the same for both models. The distribution of induced electrical potential in the VC was simulated using FEM [25]. All the simulations were carried out using COMSOL while considering the quasi-static approximation of Maxwell's equations. In this approximation, the tissues are considered to be purely resistive. This is valid for the low stimulation frequency used in TENS. However, the current volume density sources r ( ) are not excluded by this method. Poisson's equation governs the electrical potentials as shown in (1).
By setting current source Q j ( ) to zero in (2), this simplifies to Laplace's formulation in (3).
The current densities on the nerve fibers were calculated from this approximate electric potential for both models based on (2) by setting external current density J e ( ) to zero.
where, W defines the entire model, s shows tissue conductivity, V e represents induced electrical field in domains and J is current density in the media. A comparatively large non-conductive (σ=1e-10 S m −1 ) sphere was defined as external boundary and Dirichlet boundary condition V 0, e = ( ) on dW (outermost surface layer of the model)) was applied as an approximation of ground at infinity [8].

Neuron model
The target fibers (Aβ fibres), whose diameters follow a Gaussian distribution with a mean of μ D =12.5 μm and standard deviation of σ D =2 μm, were constructed based on experimental data [26] while the associated parameters were derived by interpolating experimental measurements [27]. The nerve fiber threshold levels were calculated using the McIntyre-Richardson-Grill (MRG) cable model of myelinated mammalian axons [27]. The nerve fiber compartments and their geometric positions along the nerve course were designed based on our previous study [22].
The electrical potentials were first simulated in COMSOL for both simplified and cellular models and then the obtained electrical potential field was exported into Neuron v7.4 [28] as extracellular electrical potential in the form of voltage pulses for each fiber. 100 fibers were considered and activation was considered as observing an action potential in the first and last nodes after the fifth pulse [8].

Results
The PAs versus the required stimulus current levels and the extracellular electrical potential variations across nerve trajectory for different electrode sizes are shown in figure 2.
In figure 2(a), the extracellular potential variations follow the same pattern for both models while the potential variations only introduce a shift along the nerve trajectory. As the electrode size is increased, this difference gets smaller and the electrical potential variation for the largest electrode size (E1) is the same for both models. The trends of PAs versus current levels are similar for largest electrodes for both models. Therefore, E1 has the appropriate size for comparing the two models.
The simplified model requires just slightly higher current levels to excite the same number of fibers compared with the cellular model for all electrode sizes, as shown in figure 2(b). The error between current thresholds for both models are inversely proportional to the electrode size. The error is lower for larger electrodes. For instance, to activate 50% of the nerve fibers, the required stimulus current errors between two models are 1.6% for the largest electrode and 6.8% for the smallest electrode size.
The current density variations on the nerve fiber based on E1 electrodes is shown in figure 3(a) for both models. Apart from the value at the peak, the current densities follow nearly the same trend for both models. The average current densities on the nerve trajectory for two models are, in turn, 0.48 and 0.475 A m −2 .
The features of the computation are compared in figure 3

Discussion
Bio-modelling is growingly becoming an essential step in the design and optimization of neuroprostheses. In such models, the electrical potential field is simulated in a VC and is then exported into a cable model as  Comparison of the features of computation for simplified (S) and cellular (C) models based on maximum electrode size (E1), t D and t s represent discretization and simulation time, respectively, and FE n number of finite elements. extracellular potential to predict the response of the nerve [9].
In this paper, the influence of the cellular structures on the fiber activation thresholds and computation cost was investigated using hybrid models and the results were then compared with the simplified model to reach an optimal model. The extracellular variations and PAs versus current thresholds of nerve fibers were used to quantify the difference between both models as well as electrode sizes. Our results show that both the electrical potential field and current thresholds have marginal error in the simplified model compared with the cellular model for electrodes larger than a certain size. The trend in figure 2(a) indicates that although the electrical potential along the nerve follows nearly the same trend for both models, the values for cellular models are negative at both ends of the nerve. This may be because of the mesh quality of the cellular model. The stimulus current and current density differences between the two models for E1 electrodes are within safe limits [29].
The edge effect for current density happens in highly conductive surfaces like metals. In this study, however, for this transcutaneous setting there exists a metallic contact that is in contact with a layer of a finite but not infinite conductance (part of the electrode) which is in immediate contact with the surface of the skin on its other side. Thus, the edge effect only happens in the metallic surface (as shown In figure 1(a).) but not in the electrode's surface in contact with the skin.
Although the error between the simplified and cellular models are increased when using the smaller electrode, it was observed that the slope of the PAs versus current was nearly the same for both models. The results showed that the required current levels are increased for both models when using smaller electrode sizes. This is because the smaller electrodes do not produce a field wide enough to cover the length of the nerve; thus, only the outskirts of the field will reach the nerve in that case which leads to higher required injected current to activate the nerve. The reason for the large error between the two models may be associated with localisation of the current in the cellular structure. These microscopic structures may be more sensitive to localised changes in the path of current. This may lead to a relatively lower stimulus current being required for the cellular model.
It is noted that the cell structure just below the electrode was only considered in this study due to the high computation cost involved. However, considering a larger cellular area beneath the electrodes may have an impact on the results. This should be investigated in future studies using a workstation with higher computational capabilities.
The discretization time, simulation time and the number of elements for the cellular model was considerably larger than that of the simplified model, indicating a vast difference in their respective computation cost.
Reducing the complexity of models may increase simulation efficiency by reducing simulation and discretization times. Considering the results in figures 2 and 3, this study indicates that the microscopic features have little effect on the PAs of fibers while impose a larger computation cost. On the other hand, the simplified model is computationally more efficient and has a sufficient level of accuracy at this stage of the design. Thus, it can be used to assess the effect of the neuroanatomical variations across different individuals and electrode settings with different arrangements in future investigations.

Conclusions
The quality of life can be significantly improved using neuroprostheses which are becoming a widely accepted therapeutic solution for neurological disorders. Hybrid models enable a thorough investigation of optimal parameters necessary for the design of efficient neuroprostheses. However, computational cost may limit the applicability of this approach. It was shown in this study that the effect of cellular structures of skin can be ignored and the simplified model may be used for the design and optimization of the target device.