Impact of neuroanatomical variations and electrode orientation on stimulus current in a device for migraine: a computational study

Objective. Conventional treatment methods for migraine often have side effects. One treatment involves a wearable neuromodulator targeting frontal nerves. Studies based on this technique have shown limited efficacy and the existing setting can cause pain. These may be associated with neuroanatomical variations which lead to high levels of required stimulus current. The aim of this paper is to study the effect of such variations on the activation currents of the Cefaly neuromodulator. Also, using a different electrode orientation, the possibility of reducing activation current levels to avoid painful side-effects and improve efficacy, is explored. Approach. This paper investigates the effect of neuroanatomical variations and electrode orientation on the stimulus current thresholds using a computational hybrid model involving a volume conductor and an advanced nerve model. Ten human head models are developed considering statistical variations of key neuroanatomical features, to model a representative population. Main results. By simulating the required stimulus current level in the head models, it is shown that neuroanatomical variations have a significant impact on the outcome, which is not solely a function of one specific neuroanatomical feature. The stimulus current thresholds based on the conventional Cefaly system vary from 4.4 mA to 25.1 mA across all head models. By altering the electrode orientation to align with the nerve branches, the stimulus current thresholds are substantially reduced to between 0.28 mA and 15 mA, reducing current density near pain-sensitive structures which may lead to a higher level of patient acceptance, further improving the efficacy. Significance. Computational modeling based on statistically valid neuroanatomical parameters, covering a representative adult population, offers a powerful tool for quantitative comparison of the effect of the position of stimulating electrodes which is otherwise not possible in clinical studies.


Abstract
Objective. Conventional treatment methods for migraine often have side effects. One treatment involves a wearable neuromodulator targeting frontal nerves. Studies based on this technique have shown limited efficacy and the existing setting can cause pain. These may be associated with neuroanatomical variations which lead to high levels of required stimulus current. The aim of this paper is to study the effect of such variations on the activation currents of the Cefaly neuromodulator. Also, using a different electrode orientation, the possibility of reducing activation current levels to avoid painful side-effects and improve efficacy, is explored. Approach. This paper investigates the effect of neuroanatomical variations and electrode orientation on the stimulus current thresholds using a computational hybrid model involving a volume conductor and an advanced nerve model. Ten human head models are developed considering statistical variations of key neuroanatomical features, to model a representative population. Main results. By simulating the required stimulus current level in the head models, it is shown that neuroanatomical variations have a significant impact on the outcome, which is not solely a function of one specific neuroanatomical feature. The stimulus current thresholds based on the conventional Cefaly system vary from 4.4 mA to 25.1 mA across all head models. By altering the electrode orientation to align with the nerve branches, the stimulus current thresholds are substantially reduced to between 0.28 mA and 15 mA, reducing current density near pain-sensitive structures which may lead to a higher level of patient acceptance, further improving the efficacy. Significance. Computational modeling based on statistically valid neuroanatomical parameters, covering a representative adult population, offers a powerful tool for quantitative comparison of the effect of the position of stimulating electrodes which is otherwise not possible in clinical studies.
Stimulation of occipital and vagus nerves, and sphenopalatine ganglion stimulation, delivered invasively, and stimulation of vagus and frontal nerves, and transcranial magnetic stimulation (TMS), delivered noninvasively, are examples of using neuromodulation to manage migraine. Invasive methods are used only in the most medically intractable patients due to their associated risks [8]. The induced electrical stimulation by TMS activates a mixture of neurons in the brain and may lead to an unorganized pattern of activity [9]. Moreover, the vagus nerve contains sensory and motor fibers and its stimulation may cause neck muscle contractions [10]. There is also a lack of large controlled studies to support the use of TMS and vagus nerve stimulation for the prevention of the migraine [11].
Migraine sufferers commonly report that migraine pain is centred in the regions innervated by supraorbital (SON) and supratrochlear (STN) nerves [12,13]. These sensory peripheral nerves are branches of the frontal nerve stemming from the ophthalmic division of the trigeminal nerve. STN and SON transmit frontal head pain through the trigeminal nucleus caudalis to the thalamus which is then transmitted to higher brain centres [14]. Transcutaneous frontal nerve stimulation (t-FNS) using Cefaly neurostimulator (Cefaly, CEFALY Technology, Liège, Belgium), which has been used as a basis for this study, has been shown to prevent episodic migraine [6]. Biphasic transcutaneous stimulus current pulses are applied via specially shaped flexible electrodes coated with a gel for good skin contact. The electrodes are placed across the forehead (see figures 1 and 2(a)) to stimulate the target nerve.
t-FNS using Cefaly has been tested in a double blind randomized controlled trial (n = 67) [15] in which the efficacy, tolerability and safety of the device have been studied. This study produced mixed results (50% response rate). A post-marketing survey (n = 2313) led to 53% satisfaction [16], while the most limiting factor was reported to be paraesthesia and a painful sensation [16,17]. These inconclusive results may be associated with neuroanatomical variations that lead to excessively high current levels being required [17]. Since Cefaly is patient-operated, these relatively high levels required may not be applied. The excessively high stimulus current levels may lead to the co-excitation of Aδ and C (nociceptive) fibers, resulting in painful sensation. In addition, as the electrodes are positioned near pain sensitive structures (e.g. eyes, frontal sinuses, veins and nerves), pain may be induced even at low current levels, further limiting the efficacy of the solution.
Despite all these caveats, there has been no robust investigation identifying the underlying causes of inefficacy in some cases. This is partly due to the physical limitations of studying the neuroanatomy of a statistically representative group of patients. In the computational models studied here, neuroanatomical features can be readily changed using their statistical distributions in anatomical literature. This creates a powerful tool for assessing how these variations lead to different neural responses for a given electrode setting. The computational models are composed of a volume conductor model and an advanced model of neural tissue referred to as a hybrid model. The induced electric field, due to an electrode setting, is simulated in the volume conductor model and the resulting electric potential values along the nerve are passed on to the neural model to simulate the response of the nerve.
The first step was to identify if neuroanatomical variations had any effect on the required stimulus current levels. Ten human head models were developed by varying thirteen neuroanatomical features including human head size and thicknesses of the tissue layers as well as variations in the courses of the nerve branches by considering their respective statistical distributions as reported in the literature. In each case, the required stimulus cur rent levels were simulated. The findings of this paper show that the combined neuroanatomical variations have a significant effect on the neural response for the electrode setting used in the Cefaly neuromodulator.
A potential improvement is to align the axis of electrodes with the target nerve, so that the electrical potential along the trajectory of the nerve changes polarity [18]. This can lead to lower required stimulus current levels [19]. This is achieved by rotating the existing electrode axis in Cefaly by 90° as shown in figure 1. This electrode alignment with the general nerve direction allows a large section of the nerves to be exposed to the stimulus current and the current spread is directed away from sensitive structures. The reduction in the required stimulus current by aligning the electrodes with the nerves is analyzed in this paper over a statistically representative set of head models.
The rest of the paper is organized as follows. Section 2 addresses the methods of generating the volume conductor model of a human head, associated tissue layers and myelinated nerve fibers. In section 3, the results of simulations related to the impact of neuroanatomical variations and electrode orientation on stimulation thresholds are presented. Key discussion points are presented in section 4. Concluding remarks and future directions are presented in section 5.

Methods
For all the subsequent simulations and operations, a computer with an Intel Core i7-6700 CPU at 3.4 GHz with 64 GB RAM was used.

Human head model validation
There is a high computational cost associated with an anatomically realistic human head model based on magnetic resonance imaging (MRI) [20,21] (figure 2(a)). It has been shown that the computational cost can be significantly reduced with marginal added error by simplifying the model in certain applications. To further reduce the difference between the two models, anatomical model was comprehensively assessed during image segmentation and meshing processes in this study. The geometrical model was constructed using the key neuroanatomical layers as shown in figure 2(b.1)) based on the anatomically realistic human head model. The results showed that the geometrical model can be used instead of the anatomical model to investigate the effects of the neuroanatomical variations and electrode orientation with sufficient accuracy but significantly lower computation cost. The methods of generating these models are discussed in the following subsections, and the reduction in the error between the two models is reported in detail in section 3.

Anatomical model
The anatomical three-dimensional (3D) model of the human head was obtained from a high-resolution MRI dataset. It is composed of 350 slices, each of which comprised of 480 × 480 pixels. Voxel dimensions were 0.5 × 0.5 × 0.5 mm for each of the x, y and z planes [21]. The MRI data was imported to simpleware ScanIP v2016.09 (Synopsys, Mountain View, USA) for image processing and data segmentation. The head tissue layers, mainly, skin, fat, muscle, eyeball, skull, cerebrospinal fluid (CSF), and the brain (gray and white matter), were segmented using both automatic and manual segmentation processes. Most of the layers were allocated a specific gray scale threshold range to facilitate the automatic segmentation process. However, manual segmentation was required for some tissue layers due to discontinuity (such as CSF), the presence of similar grayscale values in the neighboring regions and smoothing sharp edges that are away from the region of interest. If the manual segmentation process was not applied to overcome these problems, a tissue with high conductivity (such as CSF) would have substantially affected the current flow in the models [22]. Additionally, automatic segmentation errors were corrected using smoothing filters (recursive Gaussian, median and mean filters), by editing the morphology or filling cavities (dilate, erode, open and close functions in ScanIP software). Finally, Boolean operations were applied to remove any overlapping sections between the tissue layers.
The trajectories of the SONs and STN cannot be identified in the MRI scans due to their relatively small diameters (~1 mm) [12]. The 3D volume conductor models of these nerves were constructed from geometric shapes (e.g. cylinders) in ScanIP based on their average statistical distributions [23,24]. Multiple cylinders of different lengths and similar diameters were used to form the trajectories of the nerves as they travelled from the inferior skull layer through the subcutaneous plane. All cylinders associated with a specific nerve trajectory were unified using Boolean operations and subsequently smoothed. It was assumed that the human head is symmetrical. Therefore, the nerve models were generated only for the left side as shown in figure 2(a.2). The electrode model was generated separately and merged with the skin. The full list of segmented anatomical tissue layers and the electrode patch are shown in figure 2(a).

Geometrical model
It has been shown that in simulations focusing on the frontal side of the head, current density is not changed significantly beyond the inferior skull layer, creating a bioelectricity boundary [25]. Furthermore, in this study, the inter-electrode spacing is relatively small, leading to a relatively localized electric field. Based on simulations, the field reached 3.3% of its initial value only 1.2 cm away (within the skull layer) from the electrodes in the realistic model. Therefore, the simplified model may be used as also demonstrated in [26]. It should be noted that simplification may not be possible if electrode spacing is relatively large as in [27]. The volume conductor model of human head tissue layers and Cefaly patch electrodes were constructed from geometric shapes in COMSOL Multiphysics v5.2a (COMSOL, Ltd, Cambridge, UK) which will be discussed in section 2.3. The head model consisted of six concentric spheres to represent skin, subcutaneous tissue, muscle, skull, CSF and brain as detailed in figure 2(b). The average thicknesses of tissue layers in the anatomical model were used to construct the layers in the geometrical model. The white and gray matters were merged and modeled as a homogeneous volume representing the brain. This is a valid approximation as the electrical potential field decays considerably before reaching them. Thus, the trajectory of the nerve was only considered within the skull to subcutaneous plane level. The emerging and termination points of the nerve were chosen to be sufficiently away from the electrode such that the electric potential would tend to zero at both ends. The same nerve trajectories were used in both of models for fair comparison as shown in figure 2. It was assumed that the distance from the midline of the head to the exit points of the nerve trajectories was the same for both models. Also, the transition points of the nerves at each of the tissue layers were identical for both head models. Since the stratum corneum (SC) layer is thin compared to other tissue layers, it was modeled as a boundary condition at the outermost boundary of skin, defined in COMSOL as 'contact impedance' with its typical thickness in both models.
It was shown, in [28] that the effect of the cellular structures on the stimulus current levels can be neglected.
The geometrical (spherical) head model, shown in figure 2(b), was used in the rest of the study instead of the anatomical head model to investigate the effect of the neuroanatomical variations and electrode orientations with less computation cost and sufficient accuracy.

Neuroanatomical variations
Ten statistically relevant human head models were designed by varying the statistical distribution of key anatomical layers and nerve branches in MATLAB. Each model with neuroanatomical layers was constructed in COMSOL for simulation. These steps are detailed in the following subsections.

Anatomical layer variations
The variations of circumference of the human head was generated based on the data from 244 subjects in [28]. The age range was 17-97 years old. To automatically generate samples of this population, the available data were converted to a normal distribution based on the average and standard deviations of the head circumferences of these subjects in MATLAB v.R2018b (MathWorks, Inc., Natic M, USA). The statistical variation of the thickness of each of the 2) the relative size of the head models to the containing boundary sphere limit is also shown. tissue layers was calculated from its available mean (µ) and standard deviation (std), (see figure 3). Since the electrodes interface with the skin in the forehead area, the statistical variations of this region were identified assuming that male skin (δs) is 1 mm thicker than female [29]. There is no direct data about the subcutaneous layer thickness (δf) on the healthy human face in the literature [30] as the subcutaneous tissue is anatomically placed between skin and muscle layers. The average thickness of this tissue was calculated from the difference between the average thickness of the skin and the average depth of the frontalis muscle. It was assumed that the statistical variations of the frontalis muscle (δm) is similar to the statistical variations of the corrugator supercilii muscle (CSM) [31]. SONs and STN exit variations represent the distance away from the forehead in the z-direction (see figure 2) while SONs branching variations represent the distance away from the nerve exit point from the skull in y -direction.
Since the outermost layer of the geometrical human head was modeled as skin, the radius of this layer was calculated from the normal random distribution of the human head circumference for each model in MATLAB. To provide more statistically distributed results, the minimum difference between the circumferences of any of the head models was set to 10 mm. The remaining neuroanatomical layers were constructed by following the statistical distributions (shown in figure 3) to complete the geometrical head model. It is noted that brain was set as the innermost sphere immediately after the skull layer. Ten different human head models were generated based on these features.

Nerve variations
The SON exits the superior orbital rim and travels beneath the CSM, eventually emerging into the plane beneath the frontalis by piercing through the CSM directly or traveling beneath the CSM and emerging from the superior border of the muscle. It is generally terminated by three superficial branches as it travels in the frontalis muscle. These are superficial medial (SONs-M), intermediate (SONs-I) and lateral (SONs-L) branches. The thread-like finer nerve strands are not considered for simplification. These three branches pierce the frontalis muscle at various points and trajectories and follow different variations in the subcutaneous planes, providing sensation to the forehead [23,24]. The STN passes through the CSM either as a single nerve branch or divides into two branches when it exits from the orbital rim. Then, it pierces through the frontalis muscle to reach the subcutaneous plane [32].
It has been shown that in the variations of the exit points, the course and the branching of the STN and SONs in different planes, vary for different individuals. Thus, to explore the impact of the variations on the percentage activations (PAs) of nerve fibers for a given stimulus current, these variations should be considered. The statistical variations of the course of the associated nerve branches and their exit points from the skull have been well documented in anatomical studies [23,24]. The SONs-I was considered as a central branch.
Since the SONs branches transitioned from the skull through the subcutaneous plane, five important transition points (from i 1 to i 5 , illustrated in figure 4(a)) were selected along the trajectory of the SONs-I. Although the position of the branching along the center branch varies, it generally occurs in the frontal muscle plane. Thus, to identify the variations of the SONs-M and SONs-L, the SONs-I was chosen as the reference point. The coordinates of the branching points (m 1 , m 2 , m 3 for SONs-M and l 1 , l 2 , l 3 for SONs-L, shown in figure 4(a)) were calculated based on anatomical studies in the literature [23]. The length of the trajectory of the SONs in the skin and subcutaneous tissue layers was assumed to be 15 mm. The anatomical transitions of the two branches of the STN are generally similar while travelling from the skull layers through the subcutaneous plane. It was observed that the stimulus current level is identical for these branches based on the PAs [20]. Thus, STN was considered as a single branch in this study. The trajectory and the transition points (n 1 to n 4 ) of the nerve are shown in figure 4(a).
To generate different anatomical patterns of nerve variations, matrix A was implemented.
Here r represents the radius of the associated anatomical layer (e.g. rm: radius of muscle), α is the uniform random variation. Different mean (µ) and standard deviation (std) are statistical variations of each point along the nerve trajectory Since the geometrical head is a sphere, the radius (r, which represents x variation in Cartesian coordinate) of the different patterns of the nerve trajectory were calculated based on the first column of A. Then, the other spherical coordinates (ϕ, θ) were generated based on random distribution of each point of y and z in A. To obtain remaining Cartesian coordinates (y , z) of each point of nerve trajectory, (A.1) and (A.2) were applied as shown in appendix.
To generate realistically smooth nerve models and prevent any zigzags in the trajectory, three criteria were enforced: where i represents the point belonging to a nerve and k indicates the transition point number. Criteria (1) and (2) were applied for the first two points of the nerves for all models. If the difference between these points was negative, the rest of the points must have met the criterion in (1). However, if the difference was positive, the rest of the nerve points must have met the criterion (2). To ensure the trajectory of the nerve originated from the orbital rim and travelled through the upper section of the forehead, criterion (3) was applied.
The criterion in (1) implies that the nerve trajectory starts from the orbital rim and travels through the midline of the forehead in the z direction. On the other hand, the criterion (2) was applied to ensure the nerve travelled away from the midline of the forehead in the z direction. The incremental trajectory of the nerve in the y direction is provided by the criterion  (3). Since the nerve trajectory was expected to have a bend in all variations, the identical consecutive points (i z k+1 − i z k = 0) were not considered. These criteria were considered for SONs-I nerve points. The same criteria and conditions were implemented for all SONs branches and STN in this study.
Note that the head size and anatomical layers were generated based on their statistical distributions and the transition points of each nerve was randomly generated using an automated variable generator script in MATLAB. This process was repeated for all models in the sequence.
Once a nerve trajectory was generated, a difference of at least 1.25 mm between x points, 1.5 mm between z points and 10 mm between y points was imposed in each case. These limits were chosen while considering practical implications with regard to the computation time. Ten variations of each of SONs branches and STN were generated, leading to a total of 40 statistically relevant group of trajectories, shown in figure 4(b). The transition points of each nerve trajectory were interpolated in MATLAB to obtain a smooth and realistic nerve trajectory.
The variations of twelve anatomical features are shown in figure 3, in which the actual distribution is shown by blue traces for five standard deviations and the red dots show the values of samples of those features in the models. The methods described above in generating the head and nerve models have clearly led to a statistically diverse model set. The completed head models were named Model 1 to Model 10 (M1 to M10).

Completion of human head modelling
Once ten statistically relevant human head sizes, anatomical layers and nerve trajectories were generated, the variations for each model were exported to COMSOL to generate complete 3D volume conductors of the head models (M1 through M10). The outermost layer of the human head (skin) was constructed from a sphere based on the circumference of the human head for each model. Other anatomical layers were modelled in COMSOL based on their statistical thickness variations generated in MATLAB. Each model's specific nerve fiber trajectory was exported to COMSOL as interpolation curves. The trajectories of the nerves (3D) were generated using either sweep or Loft tools in COMSOL. Sample nerves are shown in figure 4(c). To ensure the electrode patch has appropriately shaped contact with the skin layer, an ellipsoidal geometric shape (a-axis: 94 mm, baxis: 80 mm, c-axis: 30 mm) was used to generate the electrode patch using the Boolean operation. These processes were repeated for the modelling of ten different head models.

Electrodes and their orientation
The electrodes are modeled based on the Cefaly design shown at the top in figure 4(c). The metal electrodes are covered with conductive gel (Comfort gel A hydrogel) to maintain good contact with the skin. The self-adhesive electrodes are placed across the forehead to stimulate the branches of the frontal nerve [6] with biphasic currents (see figure 1).
To assess the effect of electrode orientation on the activation current thresholds of nerve fibers, the electrode orientation was rotated by 90° for each head model as shown figure 4(c). The same size standard Cefaly electrodes (L × W = 94 mm × 30 mm) were used for the vertical electrodes for fair comparison.
To ensure appropriate placement, the vertically orientated electrodes were shifted over to the closest average nerve (STN) [24]. Note that the same thickness (Cefaly electrode thickness is 1.5 mm) and the same distance (4 mm) between source and sink were used for both electrode orientations. Since the metal layer is thin compared to other layers, it was modeled as a boundary layer at the outermost boundary of the gel layer for both electrode orientations.

Volume conductor simulations
The finite element method (FEM) was used to compute the electrical potential distribution for volume conductors. Model domains were discretized using tetrahedral finite elements for numerical solutions of partial differential equations in COMSOL. Since models are composed of geometric shapes, the domains can be meshed with finer discretization settings to obtain optimum mesh quality (without increasing computation cost significantly). Thus, electrode, skin, subcutaneous tissue and muscle were meshed using a minimum element size of 1 µm and relatively lower growth rate (1.2) and the remaining tissue layers were meshed with relatively higher minimum element sizes (e.g. 0.1-1 mm) to obtain results in a reasonable time. Meshing for the nerve model is challenging as the nerve passes through multiple layers. The mesh settings for the nerves were adjusted by different sizes and growth rates for different models. Since the outermost layer (sphere) was far from the region of interest, the discretization element size was selected to be larger (Normal tetrahedral setting) than the neuroanatomical layers. The number of elements varied approximately between seven and ten million during the discretization process, depending on the model.
As the repetition rate of the applied pulses are relatively low, the simulation results can be obtained by means of a quasi-static approximation of Maxwell's Equations in COMSOL that can be expressed by Laplace equation, as in (4). In this approximation, the tissues are considered purely resistive with no capacitive and inductive effects.
where σ is the low frequency conductivity of each of the tissues (i.e. considered with no frequency dependent elements), V FE is the electrical potential in the representative geometry.
Since zero electrical potential is defined at infinity in physics, a sphere with a large diameter (50 cm, see figure 2(b.2)) was selected around the volume conductor model to represent this boundary. The source to sink the current step was set to 1 mA.
The Dirichlet boundary condition (5) was applied to the external boundaries of the sphere. This approximates to ground at the infinity boundary condition.
where δΩ represents the outermost surface layer of the model. The conductivities of the components within the model are listed in table 1. It was shown in [33] that the anisotropic conductivity of the nerve can be ignored. Thus, the conductivity of the outermost layer (epineurium) of the nerve was used as the combined nerve conductivity. Anisotropic conductivity of a muscle was considered in this study.
Since frontalis muscle fibers are nearly vertically oriented in the longitudinal plane, the diagonal matrix of the conductivity was applied. The outer sphere (air) domain conductivity was set to a low value (1 × 10 −10 S m −1 ). Solutions were computed using an iterative conjugate gradient with algebraic multigrid preconditioner in COMSOL.

Myelinated nerve fiber model
The nerve fiber excitation was simulated using the McIntyre-Richardson-Grill (MRG) membrane mechanisms of a myelinated mammalian axon [34]. The number of compartments and their geometric positions along the nerve length were designed using a double layer cable model of myelinated fiber based on [33]. The compartments with passive properties between every two active nodes included two myelin attachment segments, two paranodal segments, and ten inter-nodal segments. Aβ fibers, whose diameters µ D followed a Gaussian distribution with a mean of 12.5 µm and standard deviation σ D of 2 µm, were modeled based on the experimental data in [35] while the associated param eters were derived by interpolating experimental measurements in [34]. Nodal length was fixed at 1 µm. The first node of Ranvier was randomly placed between 0 and Δx of the beginning of each nerve, where Δx is node to node distance for a given fiber. The compartments were inserted between every two active nodes along the nerve trajectory and the fiber model was terminated by a node based on the trajectory defined in the FEM model. 100 nerve fibers were generated for each nerve branch following the same process.
To generate a double layer cable model of 100 associated fibers with imperfect insulation of the myelin sheath, the geometric properties of all the fibers and the numbers of compartments and their positions along the nerve trajectory were exported to NEURON v7.4 [36]. The electric properties of the channels were obtained from [34] and the node dynamics were extracted from the files on Model DB with accession number 3810 [37].

FEM-NEURON models coupling
To quantify the PA of fibers, the electrical potential solutions of FEM were integrated with NEURON models. The electrical potential along the nerve was solved in COMSOL, interpolated in Matlab, and imported as extracellular potentials in NEURON to excite myelinated fibers. In NEURON, the extracellular potentials were pulsed as symmetrical biphasic current pulses of 250 µs repeated at 60 Hz as used in [6]. The amplitude of five identical consecutive current pulses was increased until activation occurred. This method is valid under quasi-static approximation. A nerve was regarded as activated when the fourth and fifth current pulses gave the same results with respect to the PAs [38]. The PAs were recorded based on the activation of the 100 fibers, where a fiber was considered activated when action potentials were observed in its first and last nodes.
As part of this study, the position of the lowest activation threshold was identified for each of the 40 nerve trajectories and for both horizontal and vertical electrode settings. To do so, a single axon model having an average diameter was formed, and the transmembrane potential was monitored at 5-node intervals while applying the extracellular potential as a stimulus pulse to the nerve model. By increasing the amplitude from zero, the position of first excitation was identified in each case.

Human head model validation
Using the Cefaly (horizontal) electrode orientation, matching anatomical and geometrical human head models were compared based on the stimulus current level leading to 50% activation of fibers. The difference in the current levels between the head models are shown in figure 5. The differences were 3% for SONs and 2% for STN. This compares to 6% and 3% respectively in [20].

Effect of neuroanatomical variations
The PAs of the fibers with respect to the required stimulus current levels for different nerve branches in the ten different head models (M1 to M10) using the horizontal Cefaly patches are shown in figure 6.
To activate all nerve fibers (PAs = 100%) in all head models, the required current range varies from 6.75 mA to 30 mA. However, it has been shown that activation of 50% of fibers is a sufficient threshold to reach the desired effects [37]. This threshold was chosen to assess the effects of the neuroanatomical variations on stimulus current levels of different nerve branches for M1-M10 as shown on each subplot. Figure 7 shows the effect of neuroanatomical variations on the amplitude of the stimulus current leading to PA = 50% of different nerve branches. These features include head circumference, different tissue thicknesses, nerve exit points and branching. The maximum and minimum current levels are 25.1 mA and 4.4 mA, respectively.
Since the current density is considerably decayed within the skull layer, the effect of the statistical variation of the skull was not studied. The statistical variation of the skull was used to provide different emerging points based on anatomical nerve variations.
Examination of figures 6 and 7 shows that there is a level of correlation between the PAs and human head size for some of the models. When human head size increases, the stimulus current thresholds also increase for, 60%, 50%, 50% and 60% of SONs-I, SONs-L, SONS-M and STN nerve models respectively. Interestingly, 30% of all nerve models have an inverse relationship. The largest head size (M10) is activated with the highest current levels for all neuroanatomical variations while the second largest head size is activated with relatively high current levels only in SONs-I nerve models.
The amplitude of the stimulus current for 50% PA of the nerve models show a level of linear relationship with the thickness of the skin. This can be seen in 50% of SONs-I and SONs-L and 60% of SONs-M and STN nerve models. However, in 20% -30% of nerve models, the required activation current thresholds are reduced as a result of increasing skin thickness.
There is a correlation between the thickness of the subcutaneous tissue and the required current thresholds for all nerve branches. An increase in the thickness of the subcutaneous tissue layer leads to relatively higher levels of stimulus current levels for 50% of SONs (SONs-, SONs-L, SONs-M) nerve models and 60% of STN nerve models. There is no variation between the remaining nerve models.
The results suggest that incremental variations in muscle thickness leads to higher levels of activation current levels for 60% of SONs-I and 50% of both SONs-M and STN. The activations of the remaining models do not change with muscle thickness variation. The muscle layer has a marginal effect (30%) on the stimulus current of SONs-L. In general, the thinnest skin and muscle layers require less stimulus current to activate all nerve variations.
The nerve exit points from the orbital rim with respect to the midline of the head are the same for SONs  branches, while the STN has different exit point variations. When the nerve shifted away from the midline of the forehead, the activation current levels increased for about 50% of all nerve models; however, three models show a negative linear relationship. The results indicate that the variations of the nerve branching ratio have an impact on the activation of the nerve fibers of the SONs branches. The branching ratio has a relative positive correlation for 50% of SONs nerve models. Such variations are not observed for 30% of all models.
In summary the minimum required current level is 4.4 mA to achieve PA = 50% for SONs fibers and the maximum current level is 25.1 mA to reach this PA for SONs-L fibers. This latter level may be beyond the range of currents which can be safely applied in this application.

Effect of electrode orientation
The results reported so far have been based on the horizontal placement of the Cefaly electrodes. The influence of the electrode alignment on the 50% activation current levels for different nerve fibers for the ten head models are shown in figure 8. The results of the horizontal electrode orientations are shown in figure 8(a) and the results of the vertical electrode orientations are shown in figure 8(b).
The current density range (not presented graphically) in the vicinity of the eyes from the surface of the skin to the back of the human head for the horizontally orientated electrodes is 140-20 mA cm −2 . This range is substantially lower (0.3-0.05 mA cm −2 ) when using the vertically orientated electrodes.
The trend in figure 8 shows that the vertical electrode orientation requires much lower current levels compared to horizontal orientation for all nerve variations based on 50% nerve fiber activation. The required current range to stimulate 50% of the fibers decreased from 4.4 mA-25.1 mA to 0.28 mA-15 mA using vertical electrode orientation. Using vertical electrode orientation, the required current levels to activate STN nerve fibers were significantly reduced and 90% were activated around 1 mA. There is also a substantial reduction in current levels for SONs nerve branches. The majority of the nerve models are activated by less than 5 mA, only one case (10%) required more than 10 mA.
The only similarity between the two arrangements is that the highest current levels are required to activate SONs-L when the models are compared.
The positions of the lowest threshold along each nerve trajectory for both horizontal and vertical electrode arrangements are shown in figure 9. The majority of nerve trajectories are activated when they exit from the skull as they travel through the muscle layers for the horizontal electrode arrangement. For the vertical electrode arrangement, the lowest thresholds are located along a portion of the nerve which is more superficial (subcutaneous plane) and close to cathode-anode separation. This concurs with the premise that this setting can lead to the participation of both electrodes in stimulation. This also indicates that the stimulation location is shifted away from sensitive structures (e.g. eyes) using the vertical electrode arrangement.

Discussion
Computational modeling can potentially facilitate the design and development of neuromodulators in various applications [25,[39][40][41][42]. In this paper a statistically representative set of models were developed and the effects of neuroanatomy and electrode orientation on the required stimulus current were studied using these models.
As the trajectory of the target nerves emerge from the skull and travel through the subcutaneous plane, various associated anatomic features may vary across different individuals [23,32]. The results indicate that while neuroanatomical variations have a substantial impact on stimulus current levels, no feature is dominant, suggesting that it is a combination of variations which contribute to these effects. Specific neuroanatomical features have a correlation with the required stimulus in some cases but this is not established and visible across all cases.
Considering the relative position of specific branches with respect to the electrodes, since SONs-L is generally further away from the two electrode orientations, this nerve branch generally requires higher levels of stimulus current when compared to other nerve branches based on simulations. The nerve fibers of the STN, on the other hand, are generally activated at a lower current threshold compared to SONs. Similarly, this may be due to the position of the exit point of the STN from the orbital rim and its proximity to the electrodes. The significance of variations in nerve position has been shown in [38].
It was noted that thicknesses of skin, subcutaneous tissue and muscle layer may affect the required stimulus current in some cases. This may be due the fact that the load seen by the electrodes changes as these thicknesses vary, and the flow of current may be affected.
For horizontally arranged electrodes, results based on neuroanatomical variations suggest that the required current levels to activate the target nerves (i.e. STN and SON), 4.4 mA-25.1 mA, were beyond the current delivery capabilities of the Cefaly device (1 mA-16 mA) in 20% of variations. This may be the cause of inefficacy in some cases. It was shown in this paper that in the majority of cases a relatively high current range (~10 mA to ~15 mA) is required to activate the nerve fibers for the horizontal arrangement. Referring to the clinical study using Cefaly, a ~ 50% successful response rate in migraine reduction has been reported [15]. Since the device is patient operated, such high current levels may lead to painful side-effects, further limiting patient compliance.
One of the key contributions of this paper is to significantly reduce the activation current thresholds by aligning the axis of the electrodes with the target nerves. The results showed that the stimulus current levels are considerably reduced and 90% of nerve models require approximately 5 mA to activate nerve fibers using vertical electrode orientation. With this improvement, the stimulus current levels for all cases are within the current capability of the existing device. It is notable that stimulating only one of the branches may be sufficient to reach the desired outcome. Considering that for the vertical electrode arrangement the STN branch and at least one branch of SONs are stimulated at current levels below 3.5 mA, significant reductions in the required stimulus current may be achieved. Such a reduction in the required stimulus current and the fact that the stimulation position will be shifted away from the eyes, may considerably improve patient compliance as undesirable side-effects-including induced pain-are less likely to be experienced.
Amongst the nerve branches, the required stimulus current for the STN was significantly reduced when using the vertical arrangement while SONs-L required the highest level of the stimulus current in 80% of the cases. As before, this may be attributed to the relative position of the electrodes and the aforementioned nerves.
It has been shown that spatially non-uniform electrical field distributions and nerve bending may lead to the activation of the nerve fibers at relatively lower thresholds. In [19], it was shown that nerve bending has a crucial impact on the required stimulus current levels due to the change in electric field distribution it introduces along the nerve, behaving as a secondary source of stimulation. The vertical orientation also introduces a spatial non-uniformity of the field.
The proposed alignment would require using two, possibly synchronized, stimulation channels. While this is not a substantial hardware change, it still adds to the complexity of the system.
A limitation of this study is the assumption that transient currents are sufficiently slow and that tissue parameters can be regarded as only resistive. Tissue conductivity values are frequency dependent and increases in frequency may result in different current distributions. Furthermore, deep branches and finer branches of the nerves were not considered which may have an impact on the current levels. However, we believe the overall conclusions of the paper would not change. Gate control theory, which constituted the premise of this study, is the most well-known of pain theories [43] based on experimental results previously reported [6]. Further studies may consider other theories of pain (e.g. pattern theory) and in tandem with experimental results, evaluate the way these theories compete in this case. Furthermore, the evaluation of such theories and the assessment of other phenomena, such as possible co-excitation of C and Aδ fibers, require the inclusion of such fibers in the model.
A clinical study of large sample size is required to assess the effect of electrode orientations on the stimulus current thresholds and pain levels to validate the computational results.

Conclusion
In this study, the multilayer FEM models of human head, transcutaneous electrical nerve stimulation (TENS) electrodes, and cable models of myelinated mammalian nerve fibers with MRG membrane dynamics were used to assess the effect of an elaborate matrix of variations of neuroanatomical structures and the influence of electrode orientations on the stimulus current thresholds in a device for migraine.
Results indicate that neuroanatomical variations have a significant impact on the stimulus current thresholds, but it is not possible to conclude if these thresholds solely depend on specific neuroanatomical variations. With horizontal axis electrodes the relatively high levels of the required stimulus currents observed for most patients are beyond the current capabilities of the existing device and exceed pain thresholds in some cases.
The proposed vertical axis electrode arrangement has multiple benefits including the reduction of the stimulus current levels and diversion of current spread from possible pain sensitive structures. This improvement can potentially improve the clinical outcome substantially if confirmed in the subsequent large sample based clinical studies.
Although the stimulus current thresholds are considerably reduced using vertical electrode orientation, the current stimulation levels still fluctuate with neuroanatomical variations and are high in some cases. The required current levels and patient discomfort (localized stimulation) can be further reduced using an optimal electrode array (using multiple sinks and sources) based on the same concept of maximizing electrical field variation in space. Here r represents the radius of the associated anatomical layer (e.g. rm: radius of muscle), α is the uniform random variation and δ the thickness of each tissue (e.g. δm: muscle thickness). This process was applied to all nerve branches.