Prediction of Force Recruitment of Neuromuscular Magnetic Stimulation From 3D Field Model of the Thigh

Neuromuscular magnetic stimulation is a promising tool in neurorehabilitation due to its deeper penetration, notably lower distress, and respectable force levels compared to surface electrical stimulation. However, this method faces great challenges from a technological perspective. The systematic design of better equipment and the incorporation into modern training setups requires better understanding of the mechanisms and predictive quantitative models of the recruited forces. This article proposes a model for simulating the force recruitment in isometric muscle stimulation of the thigh extensors based on previous theoretical and experimental findings. The model couples a 3D field model for the physics with a parametric recruitment model. This parametric recruitment model is identified with a mixed-effects design to learn the most likely model based on available experimental data with a wide range of field conditions. This approach intentionally keeps the model as mathematically simple and statistically parsimonious as possible in order to avoid over-fitting. The work demonstrates that the force recruitment particularly depends on the effective, i.e., fiber-related cross section of the muscles, and that the local median electric field threshold amounts to about 65 V/m, which agrees well with values for magnetic stimulation in the brain. The coupled model is able to accurately predict key phenomena observed so far, such as a threshold shift for different distances between coil and body, the different recruiting performance of various coils with available measurement data in the literature, and the saturation behavior with its onset amplitude. The presented recruitment model could also be readily incorporated into dynamic models for biomechanics as soon as sufficient experimental data are available for calibration.


I. Introduction
Magnetic stimulation is a method for activating neurons noninvasively through electromagnetic induction with strong and brief magnetic pulses. At present, magnetic stimulation focuses nearly exclusively on the brain [1]. Administered transcranially, magnetic stimulation can evoke direct effects, such as motor-evoked potentials [2], [3] or phosphenes [4], while certain pulse rhythms or patterns can also modulate neural circuits and shift their excitability with respect to endogenous signals [5]. However, the development of magnetic stimulation has been strongly related to the periphery; even the first successful experiments were performed on lower motor fibers and not the brain [6]. Long straight axons have been used for studying the basics of excitation since then [7]- [11].
In neuromuscular applications, magnetic stimulation can serve as an almost pain-free alternative to transcutaneous electrical stimulation [12]- [26]. Classical rehabilitation can be supported by evoking muscle contraction or performing more complex tasks, such as cycling [27]. In rehabilitation, neuromuscular stimulation serves to counteract muscle atrophy and to support relearning of movement sequences. Furthermore, orthodromic signals traveling from the periphery back to the central nervous system seem to trigger supportive neuroplastic effects [29].
Particularly the earlier neuromuscular magnetic stimulation approaches stimulated the major nerve trunks before these enter the muscle in the hope of achieving strong muscle contraction [16], [17], [20]. Researchers optimized coils to increase targeting and the recruitable force [21]. If targeted well, the evoked forces can be reasonable, but the handling is substantially more complicated than in electrical stimulation and requires experienced operators as the operator has to locate the nerve trunk and place a focal coil very accurately. Targeting and coil placement are further hampered by the contraction of the muscles and associated shifts in the anatomy. In any case, this method easily reaches rather high stimulation amplitudes for effective contraction (e.g., 70% to 100%, often on machines with already increased base power level), and still only a small portion of subjects or patients reach their maximum force levels [23], [24]. These high stimulation amplitudes can cause distress despite the better tolerability than electrical stimulation and necessarily cause extreme heating problems in the coils and pulse sources when used for training purposes [28].
Thus, recent research and clinical efforts prefer the stimulation of the intramuscular axon tree instead, where the procedure is more practical [12], [18], [19], [25]. Whereas the activation used to be weaker initially, appropriate coil designs could improve recruitment and demonstrate that better technologies can overcome such weakness and surpass electrical stimulation [13].
However, the design of novel technology for neuromuscular stimulation, including better coil geometries, requires a quantitative understanding of the recruitment. Currently, there is a major knowledge gap between the physics and the neurophysiology of neuromuscular stimulation. Consequently, it is also not clear which physical quantity is responsible for the neuromuscular activation. Thus, also the optimization of coils is rather ad-hoc. The technology for neuromuscular stimulation is therefore only improving slowly and falling behind the rapid developments in transcranial magnetic stimulation [30], [31].
For the brain, in contrast, recruitment of magnetic and electrical stimulation has been studied intensively [32]- [34] so that both experimental data sets and realistic recruitment models are available [35], [36]. Furthermore, such measurements allowed matching models with experimental data so that the dominant physical quantities could be identified [37]- [39]. Although the physical and neurophysiological conditions in the brain are obviously rather different to a peripheral muscle, the work on brain stimulation can still serve for inspiration.
Recent three-dimensional field modeling techniques and experiments provide the basis to study the activation problem quantitatively for the quadriceps femoris muscle and identify the relevant relationships between the field characteristics of the stimulation coil and the muscle activation [41], [43]. The experimental study tested a variety of coils that generate widely different field conditions in the thigh and generated further variations through different coil-leg distances with flexible spacers. The wide range of field conditions generated by the combination of both allowed to rule out that the gradient of the electric field plays a substantial role in the stimulation of the intramuscular axon tree, which pervades the muscle densely with fine branchlets, contacts each individual muscle fiber, winding with rather small curvature around them, and forms a high number of terminals.
However, all available analyses in the literature are mere correlation studies, which neither estimate muscle recruitment nor force generation. The available experimental data in the literature, on the other hand, form a sufficiently large data source with enough parameter variation to support the design of a digital twin. The combination of a realistic 3D anatomy model with a parametric recruitment model that estimates the force from anatomical (muscle and fiber anatomy) and physical (induced electric field distribution) output of the 3D model promises to close the gap between the physics and the force response. The data furthermore can serve for the calibration of free parameters to close the present gaps in the understanding.
The presented model estimates the recruitment behavior, which can be observed in isometric stimulation experiments. Appropriate mechanical descriptions are well known in the literature. Riener et al., for instance, propose a sophisticated implementation of Hill dynamics together with biomechanics and a first-order fatigue/recovery description [44]. However, the neuromuscular recruitment in such models is usually represented by a sigmoid fitting curve. This work will demonstrate that a realistic recruitment by neuromuscular magnetic stimulation can be estimated directly based on the field conditions to replace ad-hoc sigmoidal fitting curves. We will furthermore show that such a model (described in Section II) with as little as one individual parameter-individual maximum force per muscle cross-section while an optional second parameter can compensate some apparent position offsets of one specific coil, APL, in the the experimental data-allows matching data to subjects from previous experiments in the literature (see Section III). The other parameter -specifically the threshold electric field magnitude-could be considered rather constant among the subjects and might reflect the typical range of all healthy subjects.

A. Anatomy
A high-resolution model of the human thigh was prepared based on the visible-human data set of the US NIH [45]. The data include 70 mm color photographs of cryosections with 1 mm spacing in z-direction, which provide substantially higher resolution than tomography scans and allow a more detailed identification of tissues and particularly boundaries in between. Similar to other models in magnetic stimulation [46]- [48], the geometry consists of macroscopic regions with dedicated electrical properties and neglects microscopic structures, such as cell membranes. This common approach assures computability. The different classes of segmented tissue elements types include skin, fat, eleven muscles or muscle groups, the femur, blood vessels, and major extramuscular nerve branches, although the latter are not the stimulation target themselves.
The data were segmented with standard image processing methods in Matlab (The Mathworks, Natick (MA), USA). The femur and the muscles were identified by simultaneous three-channel analysis of the color data and segmented by thresholding with a manually fine-tuned tolerance band which was chosen accordingly. Furthermore, visible structures that delimit and adjoin the muscular tissue, such as tendons, supported the separation of different muscles along their boundaries and a reconstruction of their surface. The separation of the fat tissue was performed in two stages. A basic frame was obtained from thresholding. As is common in image segmentation, the threshold was determined in the corresponding histograms as a compromise between wrongly identifying other tissue types (false positive) and forming holes due to unclassified regions (false negative). Afterwards, the data set was cleaned by eroding and reconstructing the mask in order to eliminate image noise and sharp edges.
The data did not exhibit enough contrast at its boundaries for extracting the skin geometry from the images because the embedding gelatin diffused into the skin. As a remedy, this cover was generated artificially by spanning a thin layer of tissue on top of the virtual body which follows the shape of the surface (see Fig. 1 in the middle): The adipose layer was dilated with a three-dimensional Gaussian filter; thresholding that data set and subtracting all other segmented tissue types as well as still unclassified interior parts formed an approximately 2 mm surface layer.
The basis for the blood vessels was prepared by subtracting all already classified regions from the original raw data. From the remaining tissue, small not segmented isles which did not belong to any blood vessels were eliminated by a region-growing algorithm to identify such geometrically isolated subspaces and subsequent removal of unconnected spots. Remaining artifacts and noise were cancelled by three-channel thresholding. Interrupted connections of the grid formed by the vessels were reconstructed by cubic interpolation between the unintentionally separated branches. The major nerve branches (femoral and sciatic nerves as well as the tibial offsprings) were segmented manually from the imaging data.
The boundaries of the individual regions were smoothened with a three-dimensional Gaussian filter in order to suppress unrealistic artifacts and moiré effects. Remaining gaps were filled with the tissue type of the nearest neighbor. For the simulations, likewise performed in Matlab, only half of the segmented geometry, namely the right thigh, was meshed. The model is depicted in Figure 1 with its various components.

B. Coils
We modelled four different coil types, namely a standard circular coil (RND15), the racetrack coil RT-120 (MagVenture, Copenhagen, Denmark), the saddle-shaped design APL from [13], [41], and a figure-of-eight coil (MC-B70, MagVenture, Copenhagen, Denmark). The former three devices are taken from the experimental study of [41]. The experimental performance of the figure-of-eight coil for magnetic stimulation of the intramuscular nerve structures is not represented in the literature, but was added to predict its force potential using the model calibration to the other coils. Furthermore, it generates falsifiable force estimates that can be tested in later experiments to stress-test the model.
The wiring of the coils was extracted from X-ray images and modeled with each individual turn [49], [50]. A simplification of the coils to single-turn representations as common in the literature was avoided here (see Fig. 5).
The smaller coils (RND15, RT-120, and MC-B70) were placed with their cross-hairs at the very same location in the center of the proximal third of the right thigh, which is known to evoke the strongest responses in neuromuscular stimulation experiments for the quadriceps muscle group [13]. The larger APL coil covers almost the full thigh. The upper edge of the coil ends at the groin. All coils are laterally rotated by 5° in outward direction, i.e., to the right for the right thigh in the model. The positions are visualized in Fig. 5. We additionally modelled the coils with data with various spacers between coil and thigh in the previous comparative study in lifted positions by 5 mm, 10 mm, and 15 mm perpendicularly to the surface.

C. Physics
Due to the low back-action of the induced currents in the tissue, the coils were implemented as unmeshed wires, which determined the magnetic vector potential A through the Biot-Savart forward solution [51]. The segmented anatomy was meshed with hexahedra and solved with a quasi-static finite volume method (FVM) with more than 70 million volume cells. The FVM enables a high degree of stability and used an established decoupled formulation for stable simulation of eddy currents as detailed in the appendix [28], [41]. The electromagnetic induction was solved quasistatically for the sinusoidal current pulse of 5 kHz of the modeled device. The electrical characteristics were assigned according to established data in the literature [52] for the relevant frequency range as reported in Table I.

D. Force Recruitment Model
To this date, 3D models have concentrated on the physical side and studied the distribution of the electric field because the background of eddy currents is well defined and tangible using Maxwell's equations. However, for the initially requested reproduction of the experimentally observed effects, a physiologic description is essential. The presented model translates the physical data from the eddy-current simulation into experimentally accessible quantities, namely the isometric force level.
The muscle recruitment and the force generation with their particular features are to a large share sourced and caused by the physics in combination with the muscle anatomy. The physics model provides the local electric fields in the individual muscles and based on those the effective activated muscle cross section or volume as described below, which in turn feed a parametric mixed model (Fig. 2). The parametric mixed model contains all open degrees of freedom and thus calibratable components of the model, whereas the remaining parts of the recruitment follow from anatomy, physiology, and physics and are implemented as described below. The parameters were subsequently calibrated to experimental data from the literature [41].
We set up several parametric models with different numbers of parameters to compare the two plausible hypotheses for force recruitment and compensation of individuality; we identified the best-suiting model through Schwarz' Bayesian information criterion, which serves as an analytical biascompensated estimate of the Kullback-Leibler divergence to trade off high-variance and high-variability in the models and identify goodness of fit [53]. The models are summarized in Table II. The activation of the intramuscular nerve tree is believed to occur in the fine structure of the axon terminals, close to the neuromuscular junctions [41], [54]. For the physiologic approximation here, a microscopic threshold is defined for every junction, which in turn also provides the threshold of the related muscle fiber. Previous research demonstrated that the primary gradient of the electric field is not of greater relevance for explaining the force generation in neuromuscular magnetic stimulation, whereas the electric field strength magnitude of the various coils correlates with the response [41], [43]. The surface effects at the axon membrane required for an excitation are very effectively generated even in homogeneous fields due to the fine fiber structure with its small curvature radii. This phenomenon appears to reflect the situation in the cerebral cortex, where the activation is also triggered by the field strength rather than by any gradient thereof [55]- [57]. Consequently, we defined a local threshold condition at position r within each muscle, which is fulfilled as soon as the norm of the induced electric field E exceeds a certain minimum value E th , ∥E(r)∥ > E th . The free threshold parameter value (see Table II) determines the x-axis of the recruitment curve and the onset of the saturation at higher stimulation amplitudes.
The recruitable force in turn was treated as an individual characteristic. We set up two fundamental models (Models 1 vs. 2 in Table II), one volume-related, one cross-sectionarea-related as follows.
Most obvious and plausible might be estimating the force based on the activated muscle volume, which was used as one model alternative (Model 2). Due to the fibrous structure of muscles, however, the force generation is not a volume-related issue, nor does the innervation support such an approach [58]. Instead, the number of parallel myofibrils determines the peak force during the onset of a contraction in time when no fatigue effects (neither short-term nor endurance effects) are notable, whereas series muscle fibers is widely irrelevant for the tendon force [59]- [61]. This number is approximately proportional to the activated cross-section area perpendicular to the fiber pennation [62], [63].
For extracting the activated share of the physiologic cross-section area of a specific muscle which fulfills the threshold condition acts as an approximation for the force, accordingly. Each muscle of the quadriceps was handled separately in the model. The corresponding cross sections in each muscle were tilted in order to reflect the pennation axis. For the force evaluation, that cross section was taken into account which led to the highest supra-threshold area. The used individual pennation angles with respect to the femur axis were 10° for the m. rectus femoris, 8° for the m. vastus lateralis, −8° for the m. vastus intermedius, and 15° for the m. vastus medialis.
The maximum force level per physiologic muscle cross section is known to be highly individual due to the dependence on the training state, the blood circulation, microphysiology, the fiber-type composition, and the actual metabolistic conditions [63]. For both fundamental model designs, we assumed that the maximum force per area is similar for all members of the relevant extensor muscle group following earlier observations and to keep the model parsimonious and therefore assigned only one parameter for muscles of the group (f i or υ i ) [63].
Thus, we set up two parametric mixed model alternatives, each with several refinement levels to be calibrated (see Table II): A first model estimated the isometric force recruitment from the physics model based on the volume activated extensor muscles, i.e., with suprathreshold electric field. The electrical field threshold was a group parameter for all subjects, the force to volume relation, which also determines the maximum recruitable force, individual (υ i for individual i). As the curvature of the APL coil was too small for some subjects, particularly when further rubber sheet spacers were inserted, a refinement allowed for an individual shift of the APL distance (APL 0,i ). The second fundamental model used the share of the activated pennation-corrected muscle cross-section (A eff ). In its simplest form, the threshold field E th was a group parameter (i.e., assumed to be a relatively constant figure, averaged out across many axon branchlets and terminals) and the force per area f i an individual parameter. Similarly, in another refinement, the APL coil was given freedom for individual distance shifts.
The force recruitment models were coupled to the 3D physics model through the electric field and the muscle anatomy (see Fig. 2). Each model was calibrated to experimental data from the literature through mixed-effects maximum likelihood regression [41]. Regression was performed through maximization of the logarithmic likelihood of the forward model (L(F = (F ij ) ij | M E th , f i , v i , AP L 0, i ), x = (x ij ) ij )) generating the measured force responses for every sample of each subject (force F ij in response to stimulus strength x ij ) by varying the parameters (E th , f i , υ i , APL 0,i ). As the experimental data set did not contain measurements of all coils in every subject, the regression routine was designed to allow for missing data but combined regression of the entire set at once with shared and individual parameters and one overall likelihood of a model; blanks due to missing data did accordingly not contribute to the likelihood but is reflected in the sample count. We evaluated Schwarz Bayesian information to serve as a model identification criterion that accounts for the different degrees of freedom of the models, particularly in the presence of group and individual parameters.

III. Results
Whereas the forward model is relatively fast, the calibration to the data as it is based on iterative optimization of the log. likelihood was computationally more demanding and needed more than 7,000 CPU-hours on a simulation server with 24 Xeon Cores and 256 GB memory for completion.
The best fitting description used threshold electric field, area-related force, and shift of the The electric field threshold across the all data amounted to E th = (70.5 ± 21.6) V/m. However, the large spread was caused by only one outlier (see Subject S09 in Fig. 3). This outlier reached 151 V/m, which may be the result of a bad fit of the curved APL coil during the experiments used here and a large effective coil-muscle distance due to a thick subcutaneous adipose layer of the specific subject, which exceeded the one in the model, rather than a really higher local threshold of the intramuscular innervation. The median threshold field was only 65 V/m. Fig. 4 depicts the results of the calibrated model with threshold and maximum force as parameters for the four different coils with various distances from the thigh. Every coil, except the figure-of-eight device, was simulated in its initial condition and for distance values of 5mm, 10mm, and 15 mm between the coil and the skin surface in perpendicular direction to reproduce the experimental data. The figure-of-eight coil was incorporated for evaluating its unaltered performance only, since experimental values for the dependence of the distance are not available in the literature. The calibration of the threshold also determines the onset of the saturation for higher amplitudes which agrees well with the experiments and provides validation (see Fig. 4). The simulation reflects the two degrees of freedom, i.e., shift of the recruitment curve with the coil-body distance and slope variations for the different coils, which were observed in experiments [41]. This is remarkable for one reason as neither of them was represented by any parameters or modelled in any way but are purely a result of the electric field distribution generated by the different coils or their movement.
The distance of the coil from the thigh shifts the threshold in an almost linear way for the observed range. The shift of the recruitment curves in the model is smaller than in experiments, particularly for the APL coil [41]. This deviation (29% for the APL coil, 12% for RT-120 in the depicted curves) could be caused by both the experimental setup, in which spacer sheets made of rubber were used. Whereas the coils in direct contact with the thigh ideally match the surface-not least because of the flexibility of the subcutaneous adipose layer-rubber spacers are rather stiff so that the distance increases to a higher extent than the thickness of the sheets. Off-standing edges and air-filled gaps in case of the APL design are difficult to be quantified and simulated correctly. In addition, the simulated model is a standard anatomy and does not represent any individual characteristics of the experimental data which are used here.
The difference in the slopes of the various coils is relatively stable in experiments as well as in the model. The standard circular coil is nearly coincident with the racetrack device RT-120. In the simulations, the APL coil presents an approximately 2.5 times higher slope in the linear range. The figure-of-eight coil-a device which is rarely used for neuromuscular stimulation of the intramuscular axon tree due to low torque but high distress compared to other devices-shows less than half the slope of the standard round coil. The evaluation of the experimental series in [41] led to a ratio of the slopes of the APL and of the standard circular coil of 2.6.
The different performance of the coils can be visualized lucidly by marking the specific subvolume which exceeds the threshold for different points on the recruiting curve. Fig.  5 illustrates the part of the quadriceps muscle (dark) above the threshold field strength. Around the threshold, both coils exhibit just marginal activation. Whereas the increase is rather slow and locally confined for the round coil, the APL device is adjusted to the shape and the size of the quadriceps; the coil even reaches saturation within the output range of the pulse source and for reasonable power. The more homogeneous recruitment might be another, physiologic argument for using coil designs similar to the APL device. The round coil, although frequently used in neuromuscular stimulation, is not able to activate a wider range of fibers, but forms relatively local hot spots. However, these inhomogeneous strain conditions could damage relaxed fibers due to the nonphysiological strain acting on them.
Although a volume-based approach of the force estimation might be more obvious at first, the growth of the suprathreshold volume in the 3D figures for higher stimulation amplitude is much higher than the experimentally observed force increase. A direct relation of the suprathreshold volume and the force response overestimates the slope of the APL coil (factor of 3.3 instead of 2.6) and predicts a notably lower threshold for this coil compared to the other devices. The difference would correspond to a distance of almost 20 mm which was not observed in the experiments [41].

IV. Conclusion
Whereas magnetic stimulation has been too weak for daily neuromuscular applications for a long time, more efficient coils eradicated this flaw. A rather simplistic understanding in combination with heuristic trial and error were sufficient for this improvement [13] and allowed the derivation of the coupling factor to estimate the maximum efficiency level [43], [71]. This work computationally reproduced the recruiting behavior of neuromuscular magnetic stimulation for the first time. The underlying model is comparably simple, but turned out to be very capable and predicted the different recruiting of various stimulation conditions correctly. If relative torques without absolute force values are sufficient, only one parameter has to be calibrated, i.e., the local intramuscular electrical threshold field. The extraction of this single parameter from measurements entails all remaining effects, such as the experimentally found slope differences of coils and the shift of the recruitment curve for increased coil-thigh distance.
It is very likely that the model can be further simplified, for instance with respect to the anatomic resolution, in order to provide a handy tool for coil designers. In this context, it could support the urgent need of experimentalists and clinicians for adequate equipment as the performance of a design from a drawing board can be subjected to a first, but rather informative quantitative evaluation without large efforts.
Furthermore, the approach was intentionally kept as simple as possible in order to demonstrate that neuromuscular magnetic stimulation behaves macroscopically, i.e., on average notably simpler than the complexity of the microscopic neuromuscular structure might suggest. The parsimony of the model avoids a high number of parameters and factors that might facilitate over-fitting. Still, the model concurs with previous experiments and stimulation studies (see e.g., [12], [72]). The simplicity also reflects the need of coil designers in academia and industry for a flexible and usable model.
The model design and calibration became possible due to experimental data of recruitment curves under a sufficient variety of electromagnetic conditions from the literature [41]. The underlying study evaluated a series of coils and identified two degrees of freedom for changing the magnetic (and induced electric) field conditions across individuals. The essential parameters of the physiologic description-given by the average threshold of the intramuscular axon microstructure and the force per cross section-can be calibrated using measurements. All further characteristics of the recruitment curve, such as the onset of saturation, result from it. In addition, the calibration provides a value for the local electric field strength at the threshold of 65 V/m in the microstructure after excluding one outlier, which is closely related to the microanatomic conditions of the intramuscular axon tree in the neighborhood of the nerve terminals. Interestingly, the threshold value is comparable to corresponding values reported for magnetic stimulation in the brain [64]. The model is able to mimic the properties of known muscle stimulation coils quantitatively. It allows predicting the different slopes of various coils correctly and also describes the shift of the threshold with increasing spacing between coil and thigh.
In addition to serving as a tool to design and test new neuromuscular stimulation equipment in silico, the model can replace the usually predefined sigmoid functions for the recruiting curve or the behavior of another contractile element in biomechanical models [44], [65], [66]. The incorporation into such a framework furthermore provides all temporal aspects of force generation, such as force onset, fatigue and pulse repetition rate, and enables the incorporation into a control loop for functional magnetic stimulation [67], [68]. However, this step will require additional experimental data and induces further work in this field.
The quality of the underlying anatomical and geometrical model is very high and among the most-detailed in magnetic stimulation. However, even under isometric conditions, muscles change their shape during contraction. Most actual applications of neuromuscular stimulation, such as cycling, are not isometric but need force estimation for the entire motion cycle with changing knee angles and associated muscle length variations. For changing knee angles, the situation becomes even more complex though. The presented model neglects that. On the one hand, this constraint kept the model simple and was important to enable this after all very first available model to explain and predict recruitment of neuromuscular magnetic stimulation from field models and anatomy, while the predictive power of the model is sufficiently high after all. On the other hand, in a future generation, motion should be included to represent the key applications of neuromuscular stimulation in rehabilitation and muscle training. Such a step, however, will require an anatomical model including the motion and the intermittently contracted muscle or better imaging scans during such a cycle. The model closed the gap between the stimulation coil and the force generated by the muscle, which prevented systematic optimization of coils and required lengthy and costly experiments. Thus, an obvious next step for future research is the improvement of coils for neuromuscular magnetic stimulation. Due to the completeness of the model from the coil design to a muscle force prediction, initial more conventional heuristic search could with an appropriate formalism even be improved by a numerical optimization approach [73].
While the model may serve well for developing the technology further, replacing the macroscopic nature by a microscopic one including the intramuscular motor nerve tree may be another next step to improve the understanding of how and where on the local level magnetic stimulation activates nerve and in turn muscle fibers.

Appendix
The formulation used here and established previously is derived from applying the current continuity to Ampère's circuital law and introducing the electric field E through Ohm's law. After enforcing Gauß' law with Poincaré's lemma, resp. Helmholtz' theorem through magnetic vector potential with B = curl A for the magnetic flux density B and Coulomb gauge, the governing equations follow with the local electrical conductivity σ and the electrical potential ϕ. The FVM turns the differential equation into with the volume V i, j,k of hexahedral volume cell (i, j, k). Spatial discretization and firstorder finitization with secondorder error [75] of the integrals delivers Δx i, j, k Δy i, j, k Δz i, j, k f i, j, k = σ i, j, k Δy i, j, k Δz i, j, k ⋅ ϕ i + 1, j, k − ϕ i − 1, j, k 2Δx i, j, k + Δx i, j, k Δz i, j, k ⋅ ϕ i, j + 1, k − ϕ i, j − 1, k 2Δy i, j, k +Δx i, j, k Δy i, j, k ⋅ ϕ i, j, k + 1 − ϕ i, j, k − 1 2Δz i, j, k + O Δx 2 , Δy 2 , Δz 2 (4) for local variables with dimensions in the subscript. The vector potential for the excitation term f i,j,k of each cell was provided by the coil's Biot-Savart solution. The second equation was amended by natural Neumann boundary conditions for ϕ on the surface and assembled into a matrix-vector equation [28], [74]. Equations were solved by the Gauß-Seidel method. Segmented model in the resolution which was used in the simulation (in the back). In the middle image and in the front, the model is peeled in order to uncover the otherwise hidden structures. Structure of the overall model with field simulation and parametric mixed recruitment model. Regression to experimental recruitment data. Each plot displays a different subject with symbols representing the measurement and the lines the force output of the 3D field model in combination with the calibrated mixed model. The experimental data set did not include measurements of all coils in every subject. Since the differences in forces are predicted from electric field models of the coil and not influenced by any calibrated parameters, the model once calibrated can estimate the recruitment of any coil, which is shown in Fig. 4. Recruitment curves of all studied coils and conditions as predicted by the selected model. Four different distances were simulated (0 mm, 5 mm, 10 mm, 15 mm). The slope ratio between the different coils as well as the threshold-shift effect are predicted in good accordance with previous reports and not generated by any of the degrees of freedom, i.e., calibration parameters. Instead, they arise solely from the anatomy and the electromagnetic conditions of the coil geometry and therefore do not vary across subjects here. The figure-ofeight coil (MC-B70) shows a substantially lower slope, which can be quantitatively tested in the future to stress-test the model. Recruiting for a round circular coil (a) and the APL design (b) at different stimulation amplitudes (20%, 35%, and 70% of maximum stimulator output): Whereas the round circular coil shows a rather local activation pattern, the APL device leads to suprathreshold stimulation in almost the whole quadriceps muscle at 70%. The amount of the muscle volume which fulfils the threshold condition is shaded in dark color; the blue stream lines illustrate the magnetic field. Goetz