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

Abstract

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 find optimal stimulation parameters that allow for a maximum area of beneficially stimulated tissue. Finite-element simulations were performed to determine the stimulation impact of the proposed implant design and to optimize the electric field distribution resulting from sinusoidal low-frequency ($f=20\mathrm{Hz}$) electric stimulation. Optimal stimulation parameters of the electrode length $h_{\mathrm{el}}=25\mathrm{m}\mathrm{m}$ and the stimulation potential ${\phi }_{\mathrm{stim}}=0.5\mathrm{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.

1. 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].

Taking advantage of bone’s bioelectric properties, electrical stimulation for bone regeneration was attempted as early as in the 1950s with Yasuda [22] as the “pioneer” in the field of bone bioelectricity and stimulation. Since then, accelerated bone regeneration due to electrical stimulation has been widely investigated in numerous in vitro, in vivo, in silico, and clinical studies [23,24].

The naturally occurring electric fields in loaded bone contribute to a complicated signaling network involved in bone modelling and remodeling that has not yet been fully understood. However, cell experiments showed that stimulation of osteoblasts with low-frequency electromagnetic fields induces enhanced collagen synthesis [25]. Brighton et al. observed significantly enhanced DNA production following electric stimulation of bone cells [26]. Further in vitro studies could prove effects of electric stimulation on processes such as proliferation, differentiation, and migration of bone mesenchymal stem cells [27] and osteoblasts [28].

There are many ways of applying certain types of stimulating signals, for example constant, pulsed or alternating currents or fields (electric, magnetic, electromagnetic). These signals may be coupled to the tissue directly, capacitively, or inductively [29,30]. The diverse stimulation parameters applied in different studies are mainly used empirically or phenomenologically. This is because the complex physiological processes of bone remodeling are still not fully decoded. In clinical practice, low-frequency stimulation with sinusoidal electromagnetic fields between 5 and $70\mathrm{V}/\mathrm{m}$ at $f=20\mathrm{Hz}$ empirically led to improved healing results [31]. Consequently, it is assumed that electric field strengths above $70\mathrm{V}/\mathrm{m}$ might lead to harmful tissue damage due to overstimulation, whereas field strength below $5\mathrm{V}/\mathrm{m}$ would show no supportive effect on bone regeneration. The method established by Kraus and Lechner [31] has been further developed. Mittelmeier et al. proposed the “bipolar induction screw system” (BISS) [32] that became an established method and that has been applied mainly to patients with loosened hip endoprostheses or femoral head necrosis. These findings were also supported by recent in vitro studies: Cell experiments applying an AC sinusoidal signal ($f=20\mathrm{Hz}$) on human osteoblasts showed voltage-dependent enhanced differentiation of the cells [33]. Using the BISS system, electromagnetic fields in combination with additional alternating electric fields at 20 $\mathrm{Hz}$ and 700 $\mathrm{m}$$\mathrm{V}$ potential difference between the electrodes of the screw showed a positive impact on bone cell viability and differentiation [34].

As for the development of electro-stimulating implants for bone regeneration, numerical simulation is a useful tool in the design process. It allows testing of possible stimulation setups and the exclusion of unfeasible designs already at an early stage. In this context, electro-stimulating implants for a hip revision system were developed and numerically optimized [35,36,37,38]. Furthermore, electrically stimulating dental implants were designed and examined in vivo [39] and in silico [40].

However, large bone voids such as critical size defects in the mandible (lower jaw) have not been considered so far. Such defects extend over several centimeters and thus do not heal spontaneously without plastic reconstruction. Critical size defects occur mainly after partial resection of the mandible due to tumors, and are of great clinical relevance with over 3100 surgeries per year in Germany alone [41]. These defects are conventionally treated with a combination of bone replacement material and an osteosynthesis plate to keep the material in position with the help of screws. Unfortunately, this approach is frequently accompanied with complications that may make a revision surgery necessary. These complications include fracture of the plate, loosening of the fixation screws, or dehiscence [42,43,44]. To avoid such a second surgery with its additional risks for the patient, our aim is to use electrical stimulation as an approach to achieve faster bone healing and better fixation of bone and implant.

Although we already proposed preliminary numerical models of an electrically stimulated human mandible [45,46], the current work for the first time examines an electrical stimulation system for the regeneration of critical size defects in mandibular bone suitable for practical application in validation experiments. Specifically, in the current study we focus on a numerical model of a porcine, i.e., minipig mandible. This was chosen because of planned in vivo experiments with minipigs that will follow this numerical study. Compared to the human mandible, minipigs provide a similar geometry of the lower jaw. For the minipig mandible, a critical size defect is stated to be 2 cm [47].

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\mathrm{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.

2. 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^® (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.

2.1. 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 $\mathrm{m}$$\mathrm{m}$. The slice thickness was 1 $\mathrm{m}$$\mathrm{m}$ and the spacing between the slices was 0.5 $\mathrm{m}$$\mathrm{m}$. 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^® 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^® 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^® 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 $\mathrm{m}$$\mathrm{m}$ [49]. Soft tissue and skin object were uniformly remeshed with a desired edge length of 3 $\mathrm{m}$$\mathrm{m}$. 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^® 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 $\mathrm{m}$$\mathrm{m}$ in width and ca. 20 $\mathrm{m}$$\mathrm{m}$ 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 $\mathrm{m}$$\mathrm{m}$ 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_{\mathrm{el}}$ are separated by an insulator of length $55\mathrm{m}\mathrm{m}$—$2h_{\mathrm{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 $\mathrm{m}$$\mathrm{m}$ in its posterior region, here the last 10 $\mathrm{m}$$\mathrm{m}$ of the electrode are realized as a thin fixation pin of 1 $\mathrm{m}$$\mathrm{m}$ in diameter. The diameter of the remaining implant is 4 $\mathrm{m}$$\mathrm{m}$ and it is positioned approximately in the center of the defect.

The electrode length $h_{\mathrm{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).

2.2. 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 $\sigma$ and relative permittivity ${\epsilon }_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\mathrm{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. Dielectric properties (electric conductivity $\sigma$ and relative permittivity ${\epsilon }_r$) of the tissues and materials used in the simulation at $f=20\mathrm{Hz}$.

Tissue/Material	$\mathit{\sigma }$ ($\mathbf{S}$/$\mathbf{m}$)	${\mathit{\epsilon }}_{\mathit{r}}$
Cortical bone [50]	0.02	25,100
Cancellous bone [50]	0.079	4,020,200
Soft tissue (buccal mucosa) [51]	0.01	$3\times {10}^6$
Skin [50]	$2\times {10}^{-4}$	1140
Ti6Al4V [53]	$5.6\times {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 $\mathbf{J}=\sigma \mathbf{E}$ and $\mathbf{D}={\epsilon }_0{\epsilon }_r\mathbf{E}$ apply, with the current density $\mathbf{J}$, the electric field strength $\mathbf{E}$, the electric flux density $\mathbf{D}$, and the dielectric constant of vacuum ${\epsilon }_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].

2.3. 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}^{\circ }$ that would apply to a pure capacitance. The CPE is described by the equation

$$Z_{\mathrm{CPE}}=K{\left(j\frac{\omega }{{\omega }_0}\right)}^{-\beta },$$

(1)

with K being the ratio of the amplitudes of voltage and current, j the imaginary unit, $\omega$ the angular frequency $\omega =2\pi f$, and ${\omega }_0=1{\mathrm{s}}^{-1}$ a normalization frequency to account for proper units of $\mathsf{\Omega }$ for $Z_{\mathrm{CPE}}$. The parameter $\beta =0,\dots ,1$ reflects the frequency dependence of the CPE and how much it deviates from a pure capacitance ($\beta =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_{\mathrm{CPE}}$ and a charge transfer resistance $R_{\mathrm{CT}}$ [56], see Figure 2. The impedance of the tissue $Z_{\mathrm{Tissue}}$ is defined by its dielectric properties $\sigma$ and ${\epsilon }_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_{\mathrm{EDL},\mathrm{i}}$ ($i=1,2$), with

$$\begin{array}{ccc}\hfill \frac{1}{Z_{\mathrm{EDL},\mathrm{i}}}& =& \frac{1}{Z_{\mathrm{CPE},\mathrm{i}}}+\frac{1}{R_{\mathrm{CT},\mathrm{i}}}\hfill \end{array}$$

(2)

$$\begin{array}{cc}=& {\left(\frac{A_{\mathrm{meas}}}{A_{\mathrm{el},\mathrm{i}}}K{\left(j\frac{\omega }{{\omega }_0}\right)}^{-\beta }\right)}^{-1}+{\left(\frac{A_{\mathrm{meas}}}{A_{\mathrm{el},\mathrm{i}}}{R}_{\mathrm{CT}}\right)}^{-1}.\hfill \end{array}$$

(3)

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_{\mathrm{meas}}=314\mathrm{m}{\mathrm{m}}^2$ [57] and the surface of the stimulation electrode $A_{\mathrm{el},\mathrm{i}}$ used in the simulations.

The parameters ($\beta$, K, $R_{\mathrm{CT}}$) to model the electrode–tissue interface have been derived via EIS measurements of polished titanium specimens [57]. Their values from [57] are $\beta =0.95$, $K={\left(Y_0[1/F]\right)}^{-\beta }\mathsf{\Omega }={(38.95\times {10}^{-6})}^{-0.95}\mathsf{\Omega }=15453\mathsf{\Omega }$, $R_{\mathrm{CT}}=137.7\mathrm{k}\mathsf{\Omega }$. 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_{\mathrm{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.

2.4. 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\mathrm{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 $\underline{\phi }$ with $\underline{\mathbf{E}}=-\nabla \underline{\phi }$. Assuming no impressed currents, this leads to

$$\nabla ·\left(\left[j\omega {\epsilon }_0{\epsilon }_r(\mathbf{r},\omega )+\sigma (\mathbf{r},\omega )\right]\nabla \underline{\phi }\left(\mathbf{r}\right)\right)=0,$$

(4)

with the imaginary unit j, the angular frequency $\omega$, the permittivity of free space ${\epsilon }_0$, the relative permittivity ${\epsilon }_r(\mathbf{r},\omega )$, and the electric conductivity $\sigma (\mathbf{r},\omega )$. 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 ($\phi =0\mathrm{V}$) at the surface of the posterior electrode and the stimulation potential to be optimized ${\phi }_{\mathrm{stim}}\ne 0\mathrm{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

$$\underline{\phi }={\underline{\phi }}_0,{\int }_{\partial \mathsf{\Omega }}-\mathbf{n}·\mathbf{J}dS=I_0=0\mathrm{A}$$

(5)

applies a constant potential ${\phi }_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\mathrm{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 $\underline{\phi }\left(\mathbf{r}\right)$. From this quantity, the complex amplitude of the electric field strength $\underline{\mathbf{E}}=-\nabla \underline{\phi }$ can be derived. Finally, the electric field norm $|\underline{\mathbf{E}}|$ is computed that is used to rate the stimulation impact.

2.5. Optimization of the Stimulation Parameters

The stimulation parameters, i.e., the electrode length $h_{\mathrm{el}}$ and the stimulation voltage ${\phi }_{\mathrm{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:

$$5\mathrm{V}/\mathrm{m}\le |\underline{\mathbf{E}}|\le 70\mathrm{V}/\mathrm{m},$$

(6)

whereas overstimulation would correspond to $|\underline{\mathbf{E}}|>70\mathrm{V}/\mathrm{m}$ and understimulation would imply $|\underline{\mathbf{E}}|<5\mathrm{V}/\mathrm{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

$$s_{\mathrm{ben}}·\int 1-\mathrm{if}\left(\right|\underline{\mathbf{E}}|\ge 5\mathrm{V}/\mathrm{m}\&\&|\underline{\mathbf{E}}|\le 70\mathrm{V}/\mathrm{m})\mathrm{d}V,$$

(7)

$$s_{\mathrm{over}}·\int \mathrm{if}(|\underline{\mathbf{E}}|>70\mathrm{V}/\mathrm{m})\mathrm{d}V,$$

(8)

$$s_{\mathrm{under}}·\int \mathrm{if}(|\underline{\mathbf{E}}|<5\mathrm{V}/\mathrm{m})\mathrm{d}V.$$

(9)

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_{\mathrm{ben}}=s_{\mathrm{under}}=4\times {10}^{-5},s_{\mathrm{over}}=2\times {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 ${\phi }_{\mathrm{stim}}$ only, since optimizing both $h_{\mathrm{el}}$ and ${\phi }_{\mathrm{stim}}$ simultaneously was unfortunately not possible due to the rather complicated model geometry. Instead, the electrode length $h_{\mathrm{el}}$ was varied parametrically for an exemplary stimulation amplitude of 1 $\mathrm{V}$. The parameter was varied between $h_{\mathrm{el}}=10.5\mathrm{m}\mathrm{m}$ and $h_{\mathrm{el}}=27\mathrm{m}\mathrm{m}$ in increments of a few $\mathrm{m}$$\mathrm{m}$ which is admissible for the manufacturing tolerances to be expected. The “optimal” electrode length yielded $h_{\mathrm{el}}=25\mathrm{m}\mathrm{m}$ (see Section 3.2). This value was then further used during the optimization of the stimulation amplitude ${\phi }_{\mathrm{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 ${\phi }_{\mathrm{stim}}=0.2\dots 4\mathrm{V}$ with an initial value of ${\phi }_{\mathrm{stim}}=0.2\mathrm{V}$. We defined a scaling factor of $2\mathrm{V}$ to ensure that this control variable is in the order of 1, enabling the optimization algorithm to work properly.

2.6. Finite-Element Simulation

The finite-element simulations were conducted with COMSOL Multiphysics^® 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 $\underline{\phi }$. 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\times {10}^{-3}$ and a factor $\rho =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.

3. Results

3.1. Electric Field Distribution

Figure 4a shows the simulated electric field norm $|\underline{\mathbf{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_{\mathrm{el}}=25\mathrm{m}\mathrm{m}$ and ${\phi }_{\mathrm{stim}}=0.523\mathrm{V}$ (details in Section 3.2). Please note that the color legend only reaches from 5–70 $\mathrm{V}/\mathrm{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\mathrm{V}/\mathrm{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\mathrm{m}\mathrm{m}$ 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\mathrm{m}\mathrm{m}$ away from the electrode, the electric field norm has decreased below $2\mathrm{V}/\mathrm{m}$. The region with overstimulation ($|\underline{\mathbf{E}}|>70\mathrm{V}/\mathrm{m}$) is restricted to a small area directly around the electrodes.

3.2. Optimized Electrode Parameters

Figure 5 shows the volume percentage of beneficially, over-, and understimulated tissue in the defect domain in dependence on the electrode length $h_{\mathrm{el}}$ for an exemplary stimulating voltage of ${\phi }_{\mathrm{stim}}=1\mathrm{V}$.

It can be observed that the electrode length has a strong impact on the volume of stimulated tissue. Over the whole range of $h_{\mathrm{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\mathrm{V}/\mathrm{m}\le |\underline{\mathbf{E}}|\le 70\mathrm{V}/\mathrm{m}$) increases linearly, reaching a plateau at around 23 $\mathrm{m}$$\mathrm{m}$–27 $\mathrm{m}$$\mathrm{m}$. At $h_{\mathrm{el}}=25\mathrm{m}\mathrm{m}$ 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 ($|\underline{\mathbf{E}}|<5\mathrm{V}/\mathrm{m}$) generally shows an inverse behavior: It decreases from ca. 77% to ca. 35% with increasing electrode length. The volume of overstimulated tissue ($|\underline{\mathbf{E}}|>70\mathrm{V}/\mathrm{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.

Based on the optimum electrode length $h_{\mathrm{el}}=25\mathrm{m}\mathrm{m}$, the stimulation potential applied to the electrode was also optimized numerically. For the optimal stimulation amplitude, we obtain ${\phi }_{\mathrm{stim},\mathrm{opt}}=0.523\mathrm{V}$. With this parameter setting, approximately 49% of the tissue in the defect is stimulated beneficially, 49% is understimulated and 2% is overstimulated.

3.3. 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 ${\sigma }_{\mathrm{defect}}=0.079\mathrm{S}/\mathrm{m}$ (cancellous bone) to ${\sigma }_{\mathrm{defect}}=0.02\mathrm{S}/\mathrm{m}$ (cortical bone) results in a volume of beneficially stimulated tissue that is increased by ca. 21%, as can be seen in Figure 6.

However, there is also an increase in the volume of overstimulated tissue: from 2.3 to 10.6%. Hence, a reduction of the stimulation amplitude might be recommendable in the course of the physical treatment to avoid tissue damage. The current flowing through the electrodes is ca. $1.244\mathrm{m}\mathrm{A}$ in the original, optimized setting (cancellous bone). With decreasing defect conductivity, it decreases to ca. $0.71\mathrm{m}\mathrm{A}$ (cortical bone), thus leading to a reduced power consumption of the implant.

4. 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_{\mathrm{el}}=25\mathrm{m}\mathrm{m}$ and ${\phi }_{\mathrm{stim},\mathrm{opt}}=0.523\mathrm{V}$, nearly one half of the defect volume is being stimulated beneficially (Figure 6, case for ${\sigma }_{\mathrm{defect}}=0.078\mathrm{S}/\mathrm{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 $\mathrm{V}$/$\mathrm{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\mathrm{V}/\mathrm{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. ${\phi }_{\mathrm{stim},\mathrm{opt}}=0.5\mathrm{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\mathrm{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 ${\phi }_{\mathrm{stim}}=1\mathrm{V}$ and ${\phi }_{\mathrm{stim}}=2\mathrm{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_{\mathrm{el}}=25\mathrm{m}\mathrm{m}$ and ${\phi }_{\mathrm{stim},\mathrm{opt}}=0.5\mathrm{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 $\mathrm{m}$$\mathrm{A}$ 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 $\mathrm{m}$$\mathrm{A}$ [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=\frac{{\phi }_{\mathrm{stim},\mathrm{opt}}}{I}=\frac{0.523\mathrm{V}}{1.244\mathrm{m}\mathrm{A}}=420.4\mathsf{\Omega }$. 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].