Establishment of a Numerical Model to Design an Electro-Stimulating System for a Porcine Mandibular Critical Size Defect

: Electrical stimulation is a promising therapeutic approach for the regeneration of large bone defects. Innovative electrically stimulating implants for critical size defects in the lower jaw are under development and need to be optimized in silico and tested in vivo prior to application. In this context, numerical modelling and simulation are useful tools in the design process. In this study, a numerical model of an electrically stimulated minipig mandible was established to ﬁnd optimal stimulation parameters that allow for a maximum area of beneﬁcially stimulated tissue. Finite-element simulations were performed to determine the stimulation impact of the proposed implant design and to optimize the electric ﬁeld distribution resulting from sinusoidal low-frequency ( f = 20 Hz) electric stimulation. Optimal stimulation parameters of the electrode length h el = 25 mm and the stimulation potential ϕ stim = 0.5 V were determined. These parameter sets shall be applied in future in vivo validation studies. Furthermore, our results suggest that changing tissue properties during the course of the healing process might make a feedback-controlled stimulation system necessary.


Introduction
Electrical stimulation of bone has received a lot of attention in recent decades. It can be employed to improve bone healing in the case of fractures [1][2][3][4][5], non-unions [6,7], or other bone defects [8][9][10] such as those resulting after tumor resection. A further benefit could be stimulation of osseous healing capacities especially in compromised situations such as irradiated sites or patients with systemic intake of antiresorptive drugs [11,12]. The reason for the therapeutic effect of electrical stimulation is seen in the imitation of the electric fields that occur naturally in bone as a bioelectric tissue. The bioelectricity of bone has been an object of study for a long time. One example is in long bones showing the accumulation of charges when exposed to mechanical stress or strain: On the concave side negative charges accumulate and bone formation can be observed, whereas positive charges and bone resorption occur on the convex side [13]. These effects are attributed to streaming potentials [14][15][16] and piezoelectricity [17][18][19]. Also, some studies hypothesize strong interdependencies between both phenomena [20,21].
The aim of the present study is to develop and numerically examine electro-stimulating devices that directly stimulate the tissue in the defect region. In this way, no external primary coil-in contrast to the BISS method [28,32]-is necessary. This would increase the patient's comfort and ensure patient compliance. In this context, this work introduces the procedure to numerically simulate and optimize an electro-stimulating device for a critical size defect in the lower jaw of a minipig. The conducted finite-element simulation studies constitute an important first step in preparation for in vivo experiments. These experiments will be described in the future.
In our current study, we hypothesize that the proposed design of a bipolar electro-stimulating device allows for a suitable region of beneficially stimulated tissue in the lower jaw of a minipig with a critical size defect. To support this hypothesis, the stimulation impact is estimated with the help of a finite-element simulation model. In addition, we assume that the resulting numerical model enables us to determine optimal stimulation parameters such as electrode length and stimulation potential for a given stimulation frequency of f = 20 Hz to achieve a maximum beneficially stimulated area. To optimize the electro-stimulating implant with respect to best possible electric field distribution within the defect region, the finite-element method was used. Prior to simulation, a 3D model of a minipig mandible with a critical size defect, surrounding soft tissue, and an electro-stimulating implant was created. The implant parameters, i.e., the length of the electrodes and the applied voltage were optimized. In addition, we assessed the impact of varying the electric conductivity assigned to the defect and concluded important requirements for the application and further development of future electro-stimulating implants.

Materials and Methods
Here, we introduce the steps that need to be performed to establish a bioelectric numerical model of electrically stimulated biological tissue. These steps include • setting up the anatomical and technical model, i.e., segmenting computer tomographic (CT) data and computer-aided design (CAD) modelling based on the segmentation. • setting up the physical model, i.e., the anatomical and technical models are assigned their dielectric properties. • setting up the corresponding boundary value problem, i.e., the needed equations for simulating the electric field distribution. • solving the boundary value problem.
The 3D models and the simulation models created in this study are available at [48]. For the Materialise 3-matic project files (Materialise, Leuven, Belgium. http://www.materialise.com) and the COMSOL Multiphysics R (COMSOL Inc., Stockholm, Sweden. https://www.comsol.com/) models the corresponding licenses are needed. Please note that the CT data and segmentation project files are excluded from public availability.
We would like to emphasize here that all modelling steps have been undertaken with respect to the planned in vivo experiments.

Anatomical and Technical Model Generation
The geometrical model that is subject to the finite-element studies is shown in Figure 1b and was built up from the results of two modelling steps: firstly the anatomical modelling, i.e., creating the model of the defective minipig mandible and its surrounding tissue; secondly, the technical modelling, where technical components such as the electro-stimulating implant or osteosynthesis plates and screws are created. The anatomical modelling of the regarded model object, i.e., the minipig mandible, was done based on CT data of a 17-month-old male Göttingen minipig, see Figure 1a. The CT data comprised 442 slices with 512 × 512 pixels each and a pixel spacing of ca. 0.39 mm. The slice thickness was 1 mm and the spacing between the slices was 0.5 mm. By segmenting the data, i.e., assigning certain thresholds of gray values to the respective biological tissue (cortical bone, teeth, and soft tissue) a first rough anatomical model was created. For this, we used the image processing software Materialise Mimics R version 19 (Materialise, Leuven, Belgium. http://www.materialise.com).
One part of the anatomical model is the mandible (cortical bone) from which the teeth were subtracted. This was done to avoid unnecessary small details in the geometry that might lead to problems during the simulation. This is valid since the field amplitude can be easily estimated to be well below the stimulation threshold in the neighborhood of the teeth. The resulting coarse anatomical model (stereolithography (STL) file) was imported into the CAD software Materialise 3-matic R version 11 (Materialise, Leuven, Belgium. http://www.materialise.com) to be further processed. The modifications included manually filling larger holes in the geometry object that resulted from subtracting the teeth, wrapping (i.e., automatically filling smaller holes and creating a "watertight" model), smoothing of the surface, and removing other geometrical artefacts such as spikes, double or intersecting triangles, and sharp triangles. The latter modifications were performed automatically using 3-matic's "Fix Wizard". In addition, manually performed local smoothing of the surface ensured removal of further unwanted and unrealistic geometric features (little bumps and unevenness of the surface). Finally, the number of triangles describing the STL object were reduced, while preserving the mesh quality.
The soft tissue domain was created analogously: Firstly, the coarse STL file resulting after segmentation was imported into Materialise 3-matic R and the ears and the upper part of the minipig head were removed. Again, here the fields can be assessed to be negligible with respect to the stimulation threshold. Secondly, wrapping, smoothing, and repairing analogously to the bone geometry was performed. To create a skin domain, the hollow operation was used to create a shell with a uniform thickness of 2.28 mm [49]. Soft tissue and skin object were uniformly remeshed with a desired edge length of 3 mm. After importing the bone geometry, all domains, i.e., bone, soft tissue, and skin, were then combined into one final model object by creating a non-manifold assembly.
Due to the limited CT scan resolution, the cancellous bone region was not modelled based on segmentation. Instead, we modelled it as a generically shaped object (see Figure 1b) directly inside the geometry preprocessor of the simulation software COMSOL Multiphysics R version 5.3a (COMSOL; COMSOL Inc., Stockholm, Sweden. https://www.comsol.com/). Apart from that, the cancellous bone layer was neglected in bone parts far away from the so-called region of interest (ROI). The ROI is defined as the volume inside the defect domain, highlighted in red in Figure 1b.
Due to technical reasons, the geometry file of the anatomical model (consisting of the mandible and the surrounding soft tissue) and the cancellous bone firstly had to be imported into COMSOL, secondly exported as a COMSOL-internal geometry file type .mphbin, and could only then be reimported and further modified inside the simulation software. The critical size defect in the angular region of the mandible, measuring 35 mm in width and ca. 20 mm in height, was modelled within COMSOL by subtracting a cuboid from the cortical bone geometry. The defect is supposed to be filled with cancellous bone from another part of the body and will be equipped with growth factors in the in vivo experiment. Finally, so-called virtual operations were used to simplify the topological structure of the geometry, namely forming composite domains and composite faces. This allows for an easier and more regular finite-element mesh generation.
In the current study, the technical modelling of the implant design and positioning was performed directly in the geometry module of COMSOL. The technical model only comprises the electro-stimulating implant, whereas a reconstruction plate is not modelled geometrically to reduce the problem size. Instead, a boundary condition (see Section 2.4) mimics the stabilizing Ti6Al4V mesh tray around the defect domain that is supposed to hold the filling material in place. It will be 3D printed for the in vivo experiments.
The proposed cylindrical electro-stimulating implant is 55 mm long and designed in a bipolar manner, see the parametrized implant geometry in Figure 1c: Two electrodes (stimulating electrode "Electrode 1" and counter-electrode "Electrode 2") of length h el are separated by an insulator of length 55 mm-2h el . This geometry was chosen due to comparably easy manufacturing and very regular field distribution around the electrodes. Because the mandible's thickness is only around 5 mm in its posterior region, here the last 10 mm of the electrode are realized as a thin fixation pin of 1 mm in diameter. The diameter of the remaining implant is 4 mm and it is positioned approximately in the center of the defect.
The electrode length h el is a fundamental parameter influencing the electric field distribution and hence the regeneration success. Therefore, this parameter will be optimized to aim for the most beneficial electric field distribution (see Section 2.5).

Generation of the Physical Model
The physical modelling step regards assigning the physical tissue and material properties and distribution to the anatomical and technical model. In this study, it is necessary to assign the dielectric tissue and material properties, i.e., electric conductivity σ and relative permittivity ε r , to the respective model domains to simulate the electric field distribution resulting from electric stimulation of the mandible. These quantities are highly frequency-dependent and their values at the stimulation frequency of f = 20 Hz have been taken from the literature [50,51] in the case of the biological tissues (a practical online tool can be found at [52]). As described earlier, the model takes into account two types of bone: cortical bone representing the mandible and cancellous bone inside the mandible and filling the defect domain.
As for the electro-stimulating implant, the assigned biocompatible material commonly used in orthopedics and maxillofacial surgery are Ti6Al4V for the electrodes and PEEK (polyether ether ketone) for the insulator. The dielectric properties of these materials have been taken from technical data sheets [53,54]. The assigned dielectric tissue and material properties are summarized in Table 1. Table 1. Dielectric properties (electric conductivity σ and relative permittivity ε r ) of the tissues and materials used in the simulation at f = 20 Hz.

Tissue/Material σ (S/m) ε r
Cortical bone [50] 0.02 25,100 Cancellous bone [50] 0.079 4,020,200 Soft tissue (buccal mucosa) [51] 0.01 3 × 10 6 Skin [50] 2 × 10 −4 1140 Ti6Al4V [53] 5.6 × 10 5 1 PEEK [54] 10 −12 3.2 In this study, the individual tissues and materials are modelled macroscopically. Furthermore, they are simplified to be linear, isotropic, and locally homogeneous. In that case, the constitutive relations J = σE and D = ε 0 ε r E apply, with the current density J, the electric field strength E, the electric flux density D, and the dielectric constant of vacuum ε 0 . Still, it must be noted that bone is in general a highly anisotropic tissue with a hierarchic microscopic substructure [55]. Nonetheless, to employ a homogenization approach is feasible in our case. This holds especially true since only macroscopic dielectric tissue properties are available in the literature [50].

Modelling of the Electrode-Tissue Interface
In the numerical model, the electrochemical processes at the electrode-electrolyte interface between the electro-stimulating implant and the conductive tissue need to be taken into account, because part of the applied potential drops over this interface layer. These processes include the capacitive charging of the electrical double layer (non-Faradaic processes) as well as transfer of charges through the interface (Faradaic processes), as they determine the ratio of flowing current and the associated voltage drop. The electrical double layer results from the redistribution of ions in the surrounding electrolyte when interacting with the charged electrode surface. Its pseudocapacitive behavior is empirically modelled by a so-called constant phase element (CPE) [56]. This means that a constant phase difference between voltage and current exists, but generally with a larger value than the −90 • that would apply to a pure capacitance. The CPE is described by the equation with K being the ratio of the amplitudes of voltage and current, j the imaginary unit, ω the angular frequency ω = 2π f , and ω 0 = 1 s −1 a normalization frequency to account for proper units of Ω for Z CPE . The parameter β = 0, . . . , 1 reflects the frequency dependence of the CPE and how much it deviates from a pure capacitance (β = 1). The electrode-tissue interface is modelled via an equivalent circuit model that is commonly used to model simple systems: a parallel connection of the impedance of a constant phase element Z CPE and a charge transfer resistance R CT [56], see Figure 2. The impedance of the tissue Z Tissue is defined by its dielectric properties σ and ε r (ref. Section 2.2). Thus, each of the two electrode-tissue interfaces of the bipolar electrode is described by one double layer impedance Z EDL,i (i = 1, 2), with In Equation (3), the contributions to the impedance have been scaled with respect to the surface of the measurement electrode used in electrochemical impedance spectroscopy (EIS) measurements with A meas = 314 mm 2 [57] and the surface of the stimulation electrode A el,i used in the simulations.
The parameters (β, K, R CT ) to model the electrode-tissue interface have been derived via EIS measurements of polished titanium specimens [57]. Their values from [57] Please note that in EIS measurements, using the capacitance parameter Y 0 instead of the parameter K (used by Richardot and McAdams [56]) is more common.
The charge transfer resistance R CT is generally non-linear, but as we use rather small stimulation voltages we may assume its measured value to be also valid in our simulations. As the interface behavior mainly depends on the surface structure and not so much on the material itself, it can be well assumed that the values measured for titanium also give a good approximation for the titanium alloy Ti6Al4V employed in the current study.

Electro-Quasistatic Boundary Value Problem
In general, the macroscopic behavior of electromagnetic fields is described by Maxwell's equations. Assuming negligible magnetic inductive effects and propagation in the low-frequency regime ( f = 20 Hz) as is commonly valid for bioelectric phenomena, Maxwell's equations can be simplified and the so-called electro-quasistatic approximation can be applied. In this case, the time-harmonic electric field is uniquely defined by the complex scalar potential ϕ with E = −∇ϕ. Assuming no impressed currents, this leads to with the imaginary unit j, the angular frequency ω, the permittivity of free space ε 0 , the relative permittivity ε r (r, ω), and the electric conductivity σ(r, ω). Details on the derivation can be retraced in, e.g., [58][59][60].
We apply terminals with Dirichlet boundary conditions at the electrodes, with ground potential (ϕ = 0 V) at the surface of the posterior electrode and the stimulation potential to be optimized ϕ stim = 0 V at the surface of the anterior electrode ( Figure 3). Homogeneous Neumann boundary conditions are applied to the surface of the skin domain representing the insulating air. The mesh tray introduced in Section 2.1 serves as an additional stabilization for the bone but also has an impact on the resulting electric field distribution. To reduce the complexity of the model and computation time, the mesh tray has not been modelled explicitly, but only been regarded as a floating potential boundary condition at the lateral and lower boarder of the defect region ( Figure 3). The floating potential boundary condition applies a constant potential ϕ 0 on the chosen surfaces, implying that the tangential electric field is zero and the electric field is perpendicular to the boundary. The value of the resulting potential depends on the integral source current I 0 . In our study, we regard I 0 = 0 A, implying that the boundary will act as an unconnected perfect conductor. As the conductivity of the mesh tray is several orders of magnitude higher than that of the surrounding tissue, this is a valid approximation. Comparative simulations in a simplified model setup showed that the error in the electric field norm compared to a fully modelled mesh tray is less than 3% but the reduction in CPU time is 45%. Considering the mentioned boundary conditions summarized in Figure 3, Equation (4) is solved for the complex electric potential ϕ(r). From this quantity, the complex amplitude of the electric field strength E = −∇ϕ can be derived. Finally, the electric field norm |E| is computed that is used to rate the stimulation impact.

Optimization of the Stimulation Parameters
The stimulation parameters, i.e., the electrode length h el and the stimulation voltage ϕ stim need to be optimized to achieve the best possible electric stimulation in the ROI. In this context, a goal function is defined concerning the volume of beneficially stimulated tissue in the ROI. Here, beneficially stimulated means that in the considered region the following applies: whereas overstimulation would correspond to |E| > 70 V/m and understimulation would imply |E| < 5 V/m. For optimizing the stimulation parameters for a most beneficial electric field distribution, the objective is a maximum volume of beneficially stimulated tissue in the ROI. At the same time, the volume of overstimulated and understimulated tissue is to be minimized. In the Optimization module of COMSOL these goals are expressed in terms of integral objective functions Simply put, the objective functions represent the volume of beneficially stimulated (7), overstimulated (8), or understimulated (9) tissue and are only computed inside the ROI. The sum of the goal functions (7)-(9) is to be minimized to achieve the best possible stimulation impact. The goal functions are scaled with scale factors s ben = s under = 4 × 10 −5 , s over = 2 × 10 −4 so that the goal functions are in the order of one. This ensures stability of COMSOL's optimization methods. Because overstimulation is especially harmful as it might lead to tissue damage, goal function (8) is weighed with a higher factor than the other goal functions.
We consider the goal functions only in the ROI because at the current application of electro-stimulation we are not interested in an optimal fixation of the electrodes inside the residual bone. In this case, one would also consider the cortical and cancellous bone domains. Instead, we are only interested in how the electric field evolves inside a large volume of tissue to be regenerated.
The control variable is the stimulation amplitude ϕ stim only, since optimizing both h el and ϕ stim simultaneously was unfortunately not possible due to the rather complicated model geometry. Instead, the electrode length h el was varied parametrically for an exemplary stimulation amplitude of 1 V. The parameter was varied between h el = 10.5 mm and h el = 27 mm in increments of a few mm which is admissible for the manufacturing tolerances to be expected. The "optimal" electrode length yielded h el = 25 mm (see Section 3.2). This value was then further used during the optimization of the stimulation amplitude ϕ stim .
The stimulation potential applied to the electrodes has been optimized using the Optimization module of COMSOL. For this purpose, the Nelder-Mead-Simplex optimization algorithm [61,62] was used. We specified an optimality tolerance of 0.001, representing the relative accuracy in the final values of the scaled control variables. Scaling of the control variables, i.e., dividing the variables by their associated scale, ensures stability of the optimization method if the scaled control variables are in the order of unity. The derivative-free Nelder-Mead algorithm explores the design space around the current iterate by evaluating the objective function. Transformations are applied on the point of the simplex with the worst objective function value. If no further improvement of the objective function is achieved with relative increments of the scaled control variables greater or equal to the tolerance, the optimization iteration stops. In our study, the stimulation amplitude was optimized within lower and upper bounds of the control variable ϕ stim = 0.2 ... 4 V with an initial value of ϕ stim = 0.2 V. We defined a scaling factor of 2 V to ensure that this control variable is in the order of 1, enabling the optimization algorithm to work properly.

Finite-Element Simulation
The finite-element simulations were conducted with COMSOL Multiphysics R version 5.4, using the Electric Currents and Electrical Circuit interfaces of the ACDC module in order to solve the electro-quasistatic approximation of Maxwell's equations (Equation (4)). The electrode-tissue interface was modelled via an Electrical Circuit interface at each electrode. The coupling between electrical circuit and the terminal boundary condition is achieved via a so-called External-I-Terminal that applies a voltage relative to ground to the circuit node, i.e., the surface of the electrodes. For the optimization of the stimulation amplitude we used the Optimization Module of COMSOL.
As for the finite-element discretization, second-order Laplace elements were used on a tetrahedral mesh to approximate the dependent variable in the model, i.e., the electric potential ϕ. The finite-element mesh consisted of ca. 1.24 million tetrahedral elements resulting in ca. 1.68 million degrees of freedom solved for. Based on a mesh convergence study we ensured that further refining the mesh would only change the electric energy in the whole computational domain by less than 0.04% with respect to the finest mesh resolution used. Additionally, we made sure that the mesh quality, especially in the ROI, is quite high, i.e., close to 1.
The computations were performed using the frequency domain solver of COMSOL that uses the iterative solver BiCGStab (biconjugate gradient stabilized method) with a relative tolerance of tol= 1 × 10 −3 and a factor ρ = 400 in COMSOL's error estimate. This ensures that the desired tolerance would be achieved even in ill-conditioned problems. The computations were performed on a Windows workstation with 24 × 3.00 GHz CPU and 256 GB RAM. The computation time for one simulation run was about three minutes on the chosen mesh. 22 simulation runs were necessary to reach the desired optimality tolerance in the optimization study. Figure 4a shows the simulated electric field norm |E| and the electric field lines in a vertical slice of the mandible bone right through the defect region and electrodes. The field plot is shown for the optimized stimulation parameters h el = 25 mm and ϕ stim = 0.523 V (details in Section 3.2). Please note that the color legend only reaches from 5-70 V/m to emphasize the thresholds for beneficial bone stimulation. The field evolves rod-shaped around the electrodes where it achieves the desired values between 5 and 70 V/m, but also quickly diminishes with increasing distance to the electrodes, as the conductivity of the tissue is comparably high.   Figure 4b depicts the volumes of beneficially stimulated and overstimulated tissue in the defect. It can be seen that the beneficial stimulation volume reaches at least 6 mm into the tissue. As the field is not limited by the mesh tray at the sides outside of the defect, here the beneficial stimulation volume extends even further into the soft tissue. However, at a distance of roughly 10 mm away from the electrode, the electric field norm has decreased below 2 V/m. The region with overstimulation (|E| > 70 V/m) is restricted to a small area directly around the electrodes. Figure 5 shows the volume percentage of beneficially, over-, and understimulated tissue in the defect domain in dependence on the electrode length h el for an exemplary stimulating voltage of It can be observed that the electrode length has a strong impact on the volume of stimulated tissue. Over the whole range of h el the stimulated volumes vary by factors of ca. 2.2 (understimulation), 2.4 (beneficial stimulation), and 12 (overstimulation). Generally, with increasing electrode length the volume of beneficially stimulated tissue (5 V/m ≤ |E| ≤ 70 V/m) increases linearly, reaching a plateau at around 23 mm-27 mm. At h el = 25 mm the volume of beneficially stimulated tissue reaches a maximum with 54% of the defect volume being stimulated. Smaller and larger electrode lengths lead to slightly reduced volumes of beneficially stimulated tissue. The volume of understimulated tissue (|E| < 5 V/m) generally shows an inverse behavior: It decreases from ca. 77% to ca. 35% with increasing electrode length. The volume of overstimulated tissue (|E| > 70 V/m) in the considered domain is generally about one order of magnitude smaller than the volumes for beneficial and understimulation. It increases from ca. 1% to ca. 12% with increasing electrode length.

Optimized Electrode Parameters
Based on the optimum electrode length h el = 25 mm, the stimulation potential applied to the electrode was also optimized numerically. For the optimal stimulation amplitude, we obtain ϕ stim,opt = 0.523 V. With this parameter setting, approximately 49% of the tissue in the defect is stimulated beneficially, 49% is understimulated and 2% is overstimulated.

Impact of Varying Tissue Conductivity
Furthermore, the electrical conductivity of the defect domain, which is not exactly known, was varied parametrically between the extreme values of cancellous bone and cortical bone. This was done on the assumption that during the healing process the formerly soft cancellous bone develops and begins to resemble more structured and dense cortical bone [63]. Recent studies proposed the electrical properties of bone as a biomarker for bone fracture healing [64], e.g., increasing bone resistance as an indicator for bone fracture healing [65]. Hence, we may assume that the electrical conductivity in the defect domain-formerly being cancellous bone-decreases and approaches that of cortical bone during bone remodeling. The simulations show that decreasing the conductivity from σ defect = 0.079 S/m (cancellous bone) to σ defect = 0.02 S/m (cortical bone) results in a volume of beneficially stimulated tissue that is increased by ca. 21%, as can be seen in Figure 6.

Discussion
In this study, a finite-element simulation model of a minipig mandible with an electro-stimulating implant without contact to the oral cavity could be built. This model enabled us to determine optimized stimulation parameters for the electro-stimulating setup. Furthermore, we could draw conclusions on the impact of changing material properties during bone healing.
The simulation results lead to a recommendation for the stimulation parameters to be used in the in vivo experiments following this study. Applying the optimized parameters h el = 25 mm and ϕ stim,opt = 0.523 V, nearly one half of the defect volume is being stimulated beneficially ( Figure 6, case for σ defect = 0.078 S/m). By specifying strict weighting of the objective function for overstimulation, in the optimized setting only 2% of the tissue receives unfavorable electric field strengths greater than 70 V/m. It is favorable that the simulated volume of overstimulated tissue is generally much smaller than the volume of beneficially stimulated tissue. Otherwise, the related tissue damage would undo the healing impact of the electric stimulation in certain regions of the defect. However, we observed that the desired field threshold between 5 and 70 V/m is not only achieved in the bone, but also in the soft tissue (Figure 4b). The consequences of this need to be carefully studied in the in vivo experiments as the impact of such fields on soft tissue is not exactly known.
Forcing the optimization to avoid overstimulation, a relatively low optimal stimulation amplitude of ca. ϕ stim,opt = 0.5 V results. This is also favorable in terms of ensuring long lifetime of the battery being used in the planned animal experiments. Furthermore, Liboff et al. [66] observed profound electrolysis at stimulation potentials higher than 1 V. Therefore, in the planned electro-stimulating devices voltages as high as this should be avoided. In this context, comparative numerical test studies with higher stimulation amplitudes of ϕ stim = 1 V and ϕ stim = 2 V showed that the volume of beneficially stimulated tissue could be increased by 11.0% and 15.6% respectively, but these lead to a 5-fold or 11-fold increase respectively in the overstimulated volume. This brings us to the conclusion that the stimulation voltage should be monitored carefully and kept low to avoid possible tissue damage during application.
To sum up, the optimized stimulation parameters h el = 25 mm and ϕ stim,opt = 0.5 V allow for an energy-efficient beneficial electrical stimulation of roughly one half of the defect volume with only very small regions of overstimulated tissue.
The computed currents of 1.244 mA flowing between electrode and tissue may be perceptible according to [67] (reproduced in [68]); however, this value is far from the threshold for muscular reactions (5 mA [67]). Consequently, the electric currents present in the optimized stimulation setup may be acceptable for long-term stimulation. However, in the end only the in vivo experiments can reveal the actual impact.
Analyzing the consequences of changing conductivity in the defect domain allowed us to roughly estimate how the stimulation impact might change over the whole treatment time. For the in vivo experiments this might extend over 6-12 weeks. We draw the conclusion that a corresponding adjustment in the stimulation parameters is necessary to ensure proper stimulation throughout the entire therapy. Specifically, decreasing the conductivity in the defect domain ( Figure 6) revealed that with the same amplitude of the signal driving the stimulation a higher percentage of the defect region could be stimulated beneficially. What this implies is that a lower stimulation potential would be sufficient to obtain the same volume of beneficially stimulated tissue, or that the healing stimulation would reach further into the tissue, enabling a better growth of new tissue into the volume. However, substantially increased overstimulation could also be observed in this case. Therefore, the ideal stimulation unit would be feedback-controlled and able to adjust the stimulation signal during the treatment. As the flowing current decreases with the healing progress, the power consumption is reduced as well. This may allow for a longer lifetime of the battery.
For the future application in vivo-and later also in the clinic-a stimulation unit that will control and monitor the stimulation signal is under development. Regarding the design of the circuitry, the electrical impedance of the electrodes in the biological tissue is an important measure. For the optimized stimulation configuration the impedance is Z = ϕ stim,opt I = 0.523 V 1.244 mA = 420.4 Ω . Ex vivo experiments are planned to validate the numerical results. These experiments will include impedance measurements as well as measurements of the electric potential around the electrode in a porcine mandible.
Regarding the assumptions and limitations of the numerical study, there are different aspects to be noted. With respect to the dielectric tissue properties, we chose to include only the conductivity but not the permittivity. This has the advantage of simplifying the mathematical and numerical model and thus reducing computational demands. This decision is well justified since it is the conductivity which has the main impact when simulating electric stimulation of biological tissue [69]. A limitation of all such numerical studies on bioelectric effects is given by the fact that the dielectric properties of biological tissues available in the literature vary strongly among different species and also depend on the experimental conditions and the specific anatomic site [70][71][72]. Furthermore, the tissue properties depend on age [73] and health conditions of the subject: osteoporosis [74,75] will have an impact, for example. In addition, we assumed the defect to be filled purely with cancellous bone. In practice however, the cancellous bone material that is taken from another part of the body will be molded to fit into the defect, and equipped with growth factors. Hence, the cancellous bone will not be available in its original structure, thus having an unknown effect on its dielectric properties. Aside from that, we assumed each single sub-domain of the model that is representing one kind of tissue to be homogeneous and isotropic, neglecting the possible impact on the dielectric properties due to tissue heterogeneity and composition. To better capture the impact of the latter assumptions and the only vaguely known dielectric tissue properties, our future studies will include Uncertainty Quantification (UQ) methods in the numerical model. With this, the impact of the uncertain input parameters, especially the electric conductivity, on the obtained stimulation parameters could be identified. UQ methods are currently gaining in importance for this reason in numerical simulations of biological tissues [76][77][78]. In addition to the uncertainties in tissue parameters, geometrical uncertainties resulting from different jaw geometries or from manufacturing limitations of the implant could be estimated in terms of stochastic measures as well.
As the electric field decays rapidly with increasing distance to the electrodes (Figure 4a) and multiple simulations are necessary in UQ studies, future simulation models of the stimulation setup will be restricted to model only parts of the jaw geometry fully, e.g., only one half of it. This will reduce the degrees of freedom solved for and therefore reduce the computation time accordingly. Faster calculations would be especially favorable regarding UQ studies. In addition, upon success of the in vivo experiments, further numerical models including critical size defects in human mandibles will be established and analyzed.
In the current study, we neglected the complex hierarchic substructure of bone by assuming homogeneous, isotropic materials. Also, the physiological processes during bone remodeling have not been modelled. Future studies should include such microscopic details, as already quite simple studies in 2D showed notably increased electric field strengths as compared to the homogeneous case [79]. Therefore, so-called multi-scale models of electrical bone stimulation should be established. These would enable the inclusion of information from micro-scale simulations in the macro-scale model via appropriate scale-bridging techniques as presented, for instance, by Chopard et al. [80].

Funding:
This research was funded by a scholarship of the "Landesgraduiertenförderung Mecklenburg-Vorpommern" and is supported by the German Science Foundation (DFG) in the scope of the CRC 1270 "Electrically Active Implants" ELAINE.