The electromagnetic fields used in the electrical stimulation of biological samples are usually considered to be slowly varying (van Rienen et al., 2005). Thus, the magnetic field is deemed negligibly small and the electroquasistatic field equation is solved ∇ ⋅ σ ω , r + j ω ε ω , r ∇ Φ ω , r = 0 , (2) where σ is the conductivity, ɛ is the permittivity and Φ is the electric potential, which is a phasor. In general, all quantities depend on the angular frequency ω and on the respective material. However, for the system considered here, this equation can be simplified. Because the surrounding air and well plate are far less conductive than the cell culture medium, they can be accounted for by an insulating boundary condition. The electrodes are not modelled explicitly but a fixed potential is assigned to each of them as a Dirichlet boundary condition, which establishes a non-zero potential drop across the cell culture medium. Furthermore, the conductivity and permittivity of the cell culture medium can be assumed as frequency-independent up to very high frequencies far above 1 MHz (see for example (Peyman et al., 2007) where the dielectric properties of NaCl, which behaves similar as a cell culture medium (Mazzoleni et al., 1986), have been investigated). In this work also aqueous KCl solution was used, for which the same holds true (Chen et al., 2003). Given that σ ≫ ωɛ for the materials and frequency range considered here (up to 5 MHz), it suffices to solve Laplace’s equation instead of Eq. 2 ∇ 2 Φ = 0 . (3)

Note that in this case, the potential is a real-valued quantity (i.e., a phasor with a phase of zero). The FEM is a suitable method to solve Laplace’s equation on realistic geometries (Rylander et al., 2013). We used NGSolve 6.2.2102 (Schöberl, 2014), which is an open-source library for higher-order FEM built on top of the mesh generator NETGEN (Schöberl, 1997). If not stated otherwise, we used second-order Lagrange elements and a second-order geometry representation.

This work aims at the validation of the numerical simulations to ensure that a digital twin has been found. A non-invasive technique for validation is preferable because it can potentially also be used in situ. An observable that can be measured non-invasively is the current through the chamber. For a potential difference U, the current can be computed from the power dissipation P (Rylander et al., 2013) P = U I = ∫ Ω σ | ∇ Φ | 2 d Ω , (4) where Φ is the FEM solution for the electric potential. This permits to compute the resistance of the cell culture medium, which is defined as R = U/I. The current can in principle also be computed by integrating the normal component of the current density over the electrode surface. In accordance with Rylander et al. (2013), we found in preliminary numerical experiments that the current computed using the surface integration converges slower. Thus, we only report results based on Eq. 4. A quantity of interest is the electric field E, which is defined as the negative gradient of the potential Φ (i.e., is a vector field). In previous work, we found that the field is almost homogeneous in the centre of the well, where the cells are located (Budde et al., 2019), (i.e., has only one non-zero field component). Hence, we only evaluate and report the field strength (the magnitude of the field vector), which in this particular case is equal to the non-zero component of the field at the centre of the well.

For all numerical experiments, we set the voltage difference U to 1 V and the conductivity σ to 1 S m⁻¹. Thus, we will compute a reference value for the current (and resistance), which can be easily adjusted by proper scaling to match the respective experimental reality. This approach is valid because Eqs 2, 3 are linear partial differential equations. We will discuss this approach in greater detail in the Results section.

It is important to establish an estimate for the error of the numerical simulation to meaningfully interpret the results of the UQ study. Hence, we performed adaptive mesh refinement using a Zienkiewicz-Zhu error estimator Zienkiewicz and Zhu (1987) for the base geometry described in Section 2.1. From a numerical point of view, the aforementioned base geometry can be considered as the worst case. The reason for this is that this configuration features the smallest possible distance of the electrode to the dish. Thus, small elements are needed to discretise the geometry around the electrode. These elements contain a comparably large error and might need additional refinement. For other geometrical configurations, the distance of the electrode to the dish is larger and thus the numerical error is expected to be smaller. Different meshing hypotheses can be used and will eventually influence the mesh quality (Schöberl, 1997). The adaptive mesh refinement strategy on a mesh that was generated using a hypothesis to generate a very fine mesh leads to a change of the current and field strength of less than 0.01%. The deviation from benchmark results, which were obtained using the commercial FEM software COMSOL Multiphysics®V5.5, were equally small. Because we expect this numerical error to be much smaller than the possible uncertainty obtained in a UQ study, we used the above-mentioned meshing hypothesis for all computations.

To assess the accuracy and reliability of the numerical model, the uncertainties of the input parameters need to be propagated through the model. In practice, uncertain input parameters have to be identified and a probability distribution for each uncertain parameter needs to be specified. Then, the numerical model needs to be run multiple times to generate results that reflect the uncertainties of the input parameters. For this purpose, there exist two main approaches: Monte Carlo (MC) methods, which draw samples from the probability distributions and Polynomial Chaos (PC) methods, which generate a polynomial representation based on the probability distributions (a surrogate model) (Lemieux, 2009; Xiu, 2010). MC methods often require thousands of model runs to yield a reliable estimate of the model uncertainty and are thus not suitable for realistic 3D models. The PC approach with a point collocation method, which is a robust method in the context of PC UQ methods, requires usually far fewer model runs and was thus the method of choice (Tennøe et al., 2018). We used a modified version ² of the Python library Uncertainpy (Tennøe et al., 2018). The polynomial order was set to four. Statistical metrics such as the mean, variance or the Sobol indices, which express the influence of the uncertain parameters on the modelling outcome, were directly computed from the PC expansion. To estimate the 5th and 95th percentile, both 10⁴ and 10⁵ samples were drawn from the surrogate model to ensure convergence. To speed up the computations, the model runs were performed in parallel on the HAUMEA high-performance computing cluster of the University of Rostock (each computing node equipped with 2 Intel Xeon Gold 6248 CPUs with in total 40 cores and 192 GB RAM).

In Figure 1, possible error sources to be included in an UQ analysis are indicated. The assumed hypotheses for the UQ computations are summarised in Table 1. Note that we did not consider the effect of the cell culture because the focus in this work was on applications involving cells seeded in 2D culture. In 2D culture, the cells adhere to the bottom of the well in a very thin layer, which is a few micrometres thick. Such a thin layer is not expected to have any influence on the current through the chamber or the impedance, which are of interest in this work. Moreover, we considered only uniform distributions. This reflects our current knowledge of the uncertainties of the individual parameters. We would like to mention that our approach can be straightforwardly used with all probability distributions that are implemented in Chaospy including, for example, the normal distribution (Feinberg and Langtangen, 2015).

The potentiostatic DC stimulation is the stimulation method, for which the above-mentioned chamber has been designed (Mobini et al., 2016). In electrochemistry, the monitoring of the time-dependent current at fixed voltage is known as chronoamperometry (Bard and Faulkner, 2001). We performed the DC measurements using the cell neurobasal medium (details are given in Section 2.3.4) inside an incubator at 37°C. Only one electrode pair in one well was used.

The potential was applied using a laboratory power supply (Voltcraft PS 405 Pro). A digital multimeter (Voltcraft VC 404) was used to ensure a constant voltage throughout the entire experiment. We used a digital multimeter (Voltcraft VC 850) together with a Bluetooth device (Voltcraft VC 810) to record the current. The sampling interval was 1 s. We applied the current for about 10 min, then short-circuited the two electrodes until the discharging current became zero and then reversed the polarity. The current was recorded for three voltages: 1 V, 1.25 V and 1.5 V (in this order). We used different stimulation chambers for the AC and DC experiments because DC stimulation caused surface oxidation, which could have caused reduced reproducibility of AC experiments. We will discuss this in the Results section.

During all measurements inside the incubator, the temperature was recorded with a thermometer. Furthermore, the temperature inside the cell culture medium was estimated by placing a temperature sensor (DrDaq, temperature sensor DD100, PicoLog 6, Pico Technology) in an adjacent well filled with the same amount of medium.

Impedance spectra were recorded using a Gamry Reference 600+ potentiostat. The input amplitudes were set to 25 mV. In preliminary numerical experiments, we also used 50 mV and did not observe a visible difference, which indicates that the selected amplitude was chosen sufficiently small to exclude electrochemical reactions at the EEI. Unless stated otherwise, the spectra were recorded from 1 Hz to 5 MHz. The electrochemical impedance spectroscopy (EIS) measurements were carried out in a two-electrode configuration (i.e., no reference electrode was used). The EIS spectra were analysed using the open-source software ImpedanceFitter (Zimmermann and Thiele, 2021). By applying a linear Kramers-Kronig validity test (Schönleber et al., 2014), it was checked, which part of the spectrum could be successfully fitted to an equivalent circuit. Usually, only points at high frequencies greater than 1 MHz and at very low frequencies below 10 Hz had to be excluded from the analysis.

For the EIS experiments, we used both an aqueous KCl solution of known conductivity (HI7030, Hanna Instruments) and the cell culture medium for the characterisation of the chamber. The KCl solution was used at ambient conditions (25°C) and the cell culture medium as described in Section 2.3.4. The conductivity of the KCl solution at 25°C is 1.288 S m⁻¹. The conductivity of the cell proliferation medium was measured with a handheld conductivity meter (LF 325-A, Wissenschaftlich Technische Werkstätten, Weilheim, Germany) and was 1.38 ± 0.05 S m⁻¹ at 37°C.

To check if the results obtained with one electrode pair can be also used for six electrode pairs, we performed EIS measurements also using six filled wells with each 3.5 ml medium. When using six wells connected in series, the measured impedance is Z = ∑ i = 1 6 Z i ≈ 6 Z 1 , (5) where Z _i is the impedance of a single well and the Z _i are expected to be similar to the previously measured impedance of a single well Z ₁. Likewise, the impedance of six wells connected in parallel is expected to be Z = ∑ i = 1 6 1 Z i − 1 ≈ Z 1 6 . (6)

We investigated the current and voltage response to pulses with a frequency of 130 Hz, which is commonly used in DBS (Krauss et al., 2020). The pulse width was chosen as 60 µs, 200 µs, or 600 µs. Both monophasic and biphasic pulses without an interphase gap were investigated.

The voltage signal was supplied by the ISO-STIM 01D unit (NPI electronics). The current signal was measured using a 1 Ω shunt resistor and amplified using a custom-built amplifier with a gain of 10. Both signals were recorded using an oscilloscope (RTB2004, Rohde&Schwarz). Note that because of the shunt resistor, not the entire input voltage drops across the stimulation chamber. To keep the influence of the shunt resistor negligible, we chose its resistance to be much smaller than the smallest expected impedance of the stimulation chamber in the relevant frequency range.

The voltage and current responses were Fourier transformed using the fast Fourier transform (FFT) method of the NumPy package (Harris et al., 2020). The impedance was estimated by dividing the Fourier-transformed voltage signal by the current signal. This technique is also known as broadband impedance spectroscopy (Sanchez et al., 2012). The applied voltages were chosen to be 1, 1.5 and 2 V such that the current amplitude was between 1 and 10 mA. In the current-controlled mode, the current amplitude was kept fixed at 6.5 mA.

The stimulation chamber was tested with adult neural stem cells (aNSCs). aNSCs were prepared from the subventricular zone of the adult mouse brain and cultured essentially as described previously (Hermann et al., 2009; Walker and Kempermann, 2014). In brief, singularised cells were cultured as monolayer cultures on poly-l-ornithine/laminin coating in serum-free proliferation medium consisting of Neurobasal A medium (+1% glutamate, 2% B27 supplement, 1% antibiotic/antimycotic supplement (all from Thermo Fisher Scientific), 20 ng ml⁻¹ epidermal growth factor (EGF), 20 ng ml⁻¹ fibroblast growth factor 2 (FGF-2; both from Peprotech), 2 μg ml⁻¹ Heparin (Sigma)). For stimulation experiments, cells were plated on 6-well cell culture plates at a density of 38 ,000 cells/cm² in proliferation medium. After 4 days, neuro-glial differentation of aNSC was initiated by changing the medium to differentiation medium consisting of Neurobasal A (+1% glutamate, 2% B27 supplement, 1% antibiotic/antimycotic supplement (all from Thermo Fisher Scientific), 100 µm cAMP (from Sigma), and 10 ng ml⁻¹ brain derived neurotrophic factor (BDNF; from Peprotech)) for 6 days.

Here, we report only the technical details of the cell culture stimulation approach including aspects on aNSC survival. The biological effects of the stimulation represent a separate set of experiments and will be reported elsewhere.

Different stimulation protocols were assessed: short-term stimulation (current-controlled stimulation for 30 min, one hour, 2 hours; voltage-controlled stimulation for 24 h) and long-term stimulation (current-controlled stimulation for 12 h per day for 4 days (in proliferation phase) or 10 days (4 days in proliferation and 6 days in differentiation phase)). In the current-controlled mode, symmetric biphasic pulses with 6.5 mA, 130 Hz and pulse width of 60 µs were used to simulate in vivo deep brain stimulation conditions (Krauss et al., 2020) and the wells were connected in series. In the voltage-controlled mode, the same waveforms were used but with an amplitude of 1.5 V and the wells were connected in parallel.

Cell viability was tested by visual inspection in short-term stimulation conditions and by quantitative cell counting in long-term stimulation experiments using DAPI staining of cell nuclei for reliable counting. During the cell experiments, the voltages and/or currents were monitored using a RIGOL DS1000Z oscilloscope. The recordings were controlled and the data saved to a laptop using the VISA interface and a self-written Python script ³ .

We observed that the measured current did not grow linearly with the applied voltage. This would have been the expected behaviour for a circuit dominated by the ohmic resistance of the cell culture medium. Instead, the current drastically decreased with increased stimulation time. Even after about 10 minutes, the currents did not converge to a steady value. We could describe the recorded current response I by a function of the following form: I t = a t + b e − t / c , (9) where a, b, and c are positive constants and t is the time. The equation could describe a superposition of a faradaic, diffusion-limited current inversely proportional to t and nonfaradaic, capacitive current decaying with exp(−t) (Bard and Faulkner, 2001). The nonfaradaic current can, for example, be interpreted in terms of a charging of the double layer or the pseudocapacitance at the electrode surface (Merrill et al., 2005; Lasia, 2005). However, more advanced measurements would be required to unambiguously explain the observed behaviour. The behaviour of platinum during electrical stimulation is still subject of ongoing research (Hudak et al., 2017). Importantly, one cannot establish a direct relation between Eq. 9 and the cell culture medium resistance, which can be computed from the FEM solution. Thus, the current recorded at a fixed voltage cannot be used to validate numerical simulations based only on Eq. 3. Instead, local potential recordings would be required (Gundersen and Greenebaum, 1985).

Eq. 9 was fitted to the experimental data using a nonlinear least-squares method. The fitted current was in good agreement with the measured current for all voltages (Figure 4, Supplementary Figure S4, S5). Studying the two parts of Eq. 9 individually (Supplementary Figure S6) revealed that before reversing the polarity, the faradaic and nonfaradaic currents are on the same order of magnitude. The nonfaradaic current is almost constant at all times. During the measurement period, the current did not converge to this constant value. After reversing the polarity, the influence of the faradaic current is considerably smaller. It even seems as if the current after polarity reversal was a continuation of the current before polarity reversal. Because we short-circuited the electrodes and thus there should be no residual charge stored in the system before reversing the currents this observation is surprising. The result suggests an electrochemical memory of the system. This could mean that both electrodes are continuously changed during the experiment and that the state of the electrodes is not reversed when changing polarity. There might also be other reasons for the irreversibility; for example, depletion of reactive species around the electrodes (Boehler et al., 2020). Then, the composition of the medium around the electrodes could have changed and a diffusion layer, which we do not include in our simulation model, could be present. When using the stimulation chamber in cell experiments, we observed a change of the colour of the anode, most likely showing oxidation (PtO₂). Thus, we cleaned the electrodes electrochemically after each DC stimulation application by applying a higher voltage of U = 5 V for 5 min in NaCl. For the cell experiments, this ensured replicable stimulation currents.

The equivalent circuit shown in Figure 2 distinguishes between the two EEIs. However, the EEIs are in practice indistinguishable unless a reference electrode is used. Hence, the EEIs are often described by one circuit comprising a constant-phase element (CPE) in parallel with a charge-transfer resistance (Richardot and McAdams, 2002). We found that this equivalent circuit did not describe the EIS spectra well. Instead, we used a circuit that had been developed to describe platinum surface oxidation (Supplementary Figure S7); more explanations on the involved elements are given in the Supplementary Material, Supplementary Figures S8, S9) (Ragoisha et al., 2010). Then, the fitted impedance values deviated usually less than 1% from the experimental data (an example is shown in Figure 5). For some configurations, an additional lead inductance improved the fit results. This was particularly the case for the 6-well configuration. The ohmic resistance R _medium can be directly compared to the numerical simulations (Table 2). All measured values lie within the prediction intervals of the numerical simulation, which validates our model. Using different volumes of the KCl solution showed that the liquid could be modelled entirely as a resistor: the measured imaginary part did not depend on the volume, which would have been expected if the imaginary part would not be exclusively due to the EEI (Supplementary Figure S10). The cutoff frequency where the impedance changes from capacitive to resistive behaviour can be estimated to lie between 1 and 10 kHz. We will later show the impact of this quantity on the current and voltage transients. In sum, these results show that the numerical simulations can reliably predict the ohmic resistance of the culture medium while the EEI properties can only be inferred from EIS measurements.

It is known that the EEI impedance behaves nonlinearly with increasing voltage amplitude at low frequencies (i.e., less than 1 kHz) (Moussavi et al., 1994; Richardot and McAdams, 2002). Thus, we checked the impedance at the fundamental frequency (130 Hz) for increasing voltage amplitudes. Indeed, we could observe nonlinear behaviour (Supplementary Figure S11). The impedance did not change notably at amplitudes lower than 250 mV. Hence, this voltage amplitude can be used as an estimate for the limit of linearity.

Rectangular waves can be described in the frequency domain by Fourier series (see also Supplementary Material). This reveals that the frequencies used in therapeutic applications such as DBS also contain high frequencies (Gimsa et al., 2005). To obtain the frequency-domain representation of the signals, there exist two popular approaches: fast Fourier transform (FFT) of time-domain signals (Butson and McIntyre, 2005) or the use of the analytically available expressions for the Fourier series (Butenko et al., 2020). The main difference in the frequency spectra of the waveforms considered in this work is that the biphasic pulse has its main contribution at higher frequencies than the monophasic pulse (Supplementary Figure S12–14). The amplitudes of the individual frequency components of the different waveforms did not exceed the aforementioned limit of linearity for an overall pulse amplitude of 1 V. Amplitudes greater than or equal to 2 V would lead to frequency amplitudes greater than 250 mV and thus potentially non-linear responses at low frequencies.

We used the FFT approach to estimate EIS spectra from time-domain data. The impedance is given by Z ω = U ω I ω (10) with ω the angular frequency, Z the impedance, U the potential and I the current in the frequency domain. The impedance computed by the FFT is known at equally spaced frequencies. The frequency resolution (i.e., the frequency spacing) depends on the length of the time signal. We found that recording about 10 periods, which corresponds to a frequency resolution of about 10 Hz, yielded a sufficient resolution. We used a truncation method to reduce noise in the FFT spectra: only data points with a current amplitude that is at least 10% of the maximum current amplitude were considered. This approach has also been proven to be effective for numerical simulations (Butenko et al., 2019).

To construct time-domain signals under consideration of the measured EIS spectra, we used the analytical Fourier series approach. Note that charge-balancing signals such as symmetric biphasic, biphasic with delay etc. can be straightforwardly computed as the superposition of time-shifted monophasic rectangular waves.

In preliminary experiments, we found that the current-controlled monophasic waveform with KCl showed a problem of a large DC voltage offset (Supplementary Figure S16). Moreover, the EIS spectra changed after using such a waveform, which indicated a possible change of the electrodes as also mentioned in Section 3.2.1 (data not shown). Hence, we did not further consider current-controlled monophasic waveforms because potentially harmful electrochemical reactions cannot be excluded. Simple reasoning for this DC offset, which has also been reported elsewhere (Paap et al., 2021), can be given based on the results presented in Section 3.2.1. The employed monophasic pulse with an amplitude of 6.5 mA, a pulse width of 60 µs, and a frequency of 130 Hz has a DC component of 50.7 µA. We found that a DC voltage of about 1 V caused a current of only about 10 µA decaying with time (Figure 4). This explains why the current-controlled monophasic waveform required a voltage DC offset greater than 1 V.

The impedance data obtained using the FFT algorithm (Figure 6) could be well explained using the EIS results from Section 3.2.2. Nevertheless, the impedance deviated slightly from the impedance measured by EIS. Thus, we fitted the impedance again to update the parameter values of the impedance model (Supplementary Figure S7). Considering a change in all variables turned out to be an inappropriate approach. Some waveforms contain information only in a limited frequency range (Supplementary Figure S12-14). In this case, not all parameters of the impedance model could be unambiguously determined. Some fit parameters were linearly correlated. In consequence, this caused a wrong estimate of the ohmic resistance. Instead, we found that permitting changes in 1) inductance, 2) ohmic resistance and 3) double-layer capacitance sufficed to accurately describe the measured data. The increased lead inductance (evidenced by the positive phase of the measured data in Figure 6) is most likely caused by the long and unshielded wires connecting the stimulator and the stimulation chamber. Because the conductivity of the medium depends on the temperature, the deviation of the expected and observed ohmic resistance can be, for example, explained by the uncertainty of the incubator’s temperature control. The change in the double-layer capacitance was usually only a few per cent. It could be caused by an electrochemical reaction at the electrode surface (Ragoisha et al., 2010), which occurs due to the applied electrical stimulation. With this result, we established a means to update the model (Supplementary Figure S7) based on evolving data, as required for a digital twin (Wright and Davidson, 2020). Furthermore, it permits to identify (undesired) changes in the stimulation system with increasing stimulation time.

For the (re-)construction of the stimulation signals, two modes have to be considered: voltage- and current-controlled stimulation. Either the voltage or the current pulse is controlled to be a rectangular pulse. Having the parameter values to compute the impedance at hand, we estimated both signals in the frequency domain using their Fourier series representation and Eq. 10. We observed a very good agreement between theory and experiment (Figures 7, 8). Particularly, the voltage and current transients in voltage-controlled mode could be predicted by the fitted parameter values of the impedance model (Figure 7 and Supplementary Figure S17). The measured currents show the influence of the cutoff frequency, which is deemed to be one important characteristic of a stimulation electrode (Boehler et al., 2020). The signals with dominant contributions at frequencies greater than the cutoff frequency (Figure 7A and Supplementary Figure S18A) yield a more rectangular current than the signals with a longer pulse width and thus more dominant low-frequency components.

Notable deviations between the prediction and the recorded data were only observed in the current-controlled mode and when the six wells were connected in parallel. The voltage in the current-controlled biphasic stimulation set-up revealed a DC offset (Figure 8). The offset could be reduced by manually tuning the stimulator before each experiment but could still amount to about 100 mV due to limited tuning accuracy. The biphasic signal does not comprise a DC component and thus we suppose that the DC offset stems from a coupling capacitor in the stimulator (van Dongen and Serdijn, 2016). Other authors have also reported a similar observation, which was termed DC contamination (Neudorfer et al., 2021). Due to the small magnitude of the DC offset, it does not significantly affect the current through the medium and is thus not identifiable without monitoring the voltage. Furthermore, without comparison to the digital twin prediction, it could possibly be overlooked. Nevertheless, the offset might still cause continuous electrochemical reactions.

While the current-controlled monophasic pulses revealed a very large DC voltage offset (Supplementary Figure S16), we observed a negative DC current offset for voltage-controlled monophasic pulses. This was most evident when using KCl solution (Supplementary Figure S17) instead of the medium (Supplementary Figure S18). In contrast to the current-controlled mode, the DC current offset was in good agreement with the theoretical prediction (Supplementary Figure S17, 18). The impedance predicted by the model (Supplementary Figure S7) tends to infinity when the frequency tends to zero (i.e., to the DC limit). Then, the (positive) DC current component is blocked by an infinitely large DC impedance (see Eq. 10). Because the DC current component does not contribute to the current signal, a significant (negative) DC offset can be observed. In terms of electrochemistry, the DC offset indicates an infinitely large charge transfer resistance, which suggests that no significant faradaic electrochemical reactions occur (Richardot and McAdams, 2002). Due to the blocking effect, the cells would be exposed to a small field in the negative direction even when no signal is actively applied (i.e., when the input voltage is zero). For this reason, current-controlled pulses are often preferred over voltage-controlled pulses because the applied stimulation field strength is proportional to the current density but not to the applied voltage. The aforementioned DC voltage offset that appears for current-controlled monophasic pulses can then be removed using charge-balancing approaches (Paap et al., 2021), which we did not cover in this work.

In the case of the parallel connection, it turned out that the waveform slightly deviated from the expected waveform (Supplementary Figure S19). The impedance of the 6-well system connected in parallel is only about 30 Ω. Hence, the total current through the system becomes large (about 60 mA) and might negatively affect the performance of the stimulator (Tandon et al., 2009). This result highlights the importance of a digital twin for the performance assurance of the electrical stimulation device. Still, the agreement between prediction and the measured current was good (Supplementary Figure S19B). By integrating the shunt resistor (1 Ω) into the equivalent circuit model (Supplementary Figure S7) and repeating the analysis, we could study its influence. At this point, the shunt resistor did not significantly change the results, but we will discuss later a case, where the shunt resistor has to be modelled explicitly.

When connecting the six wells in series, we did not observe similar behaviour as for the parallel connection (data not shown). This indicates that the observed deviations are indeed explained by the small load impedance of the parallel connection. In general, the good agreement between predicted and measured voltage and current transients for the parallel and series connection is highly important because it suggests that there is no significant difference between the individual electrode pairs with respect to their electrochemical behaviour.

Validating the FEM simulations enabled us to establish a connection between macroscopic quantities (voltage, current) and local quantities (potential, field strength). This permits estimating the field strengths to which the cells are exposed from the current transients. For that, the voltage drop across the medium U (i.e., the boundary condition of the simulation) is computed by multiplying the measured current I and the computed resistance R (known from Eq. 8 for a known conductivity σ). We use the UQ bounds for the resistance R (Table 2) to obtain error bounds for the voltage drop U. For each U, the prediction interval for the field strength is known (see Section 3.1). Because this approach requires knowledge of the current I and the conductivity σ, it is termed current-conductivity method.

There is a second way to estimate the field strength inside the cell culture medium through the estimation of the voltage drop across the cell culture medium U _medium. This approach requires exact knowledge of the impedance of the medium Z _medium to apply the voltage divider formula U medium ω = Z medium ω Z ω U in ω . (11)

The total impedance Z is known from the fitted EIS spectra and/or fits to the Fourier-transformed voltage/current transients. This approach, which we refer to as the voltage-divider approach, has one advantage over the previously presented current-conductivity approach: the error of the conductivity does not need to be considered and thus the field estimates are more accurate (Figure 9).

The same field estimate procedure can be applied for wells connected in series or parallel. However, then the error estimates are less reliable because the computation is done using the approximation that all electrode pairs share the same impedance (see Eqs 5, 6). We attempted to account for this by a worst-case assumption: the lower bound for the impedance is computed under the assumption that all wells have the minimal impedance computed for one well. Vice versa, the maximum impedance of one well was used to estimate the upper bound for the impedance. For the voltage-divider approach, we did not include an additional error estimate and used the error bounds for the field strength determined for one well. This approach probably underestimates the error.

The presented approach for estimating the electric field strength based on an equivalent circuit model is usually referred to as lumped-element approach. Other options are distributed impedance models (Cantrell et al., 2008; Howell et al., 2014), which integrate the EEI impedance as a Robin boundary condition into the FEM model (more details are given in the Supplementary Material). Thus, they require simulation runs for each EEI impedance. Furthermore, they suffer from the problem that the EEI impedance of the individual electrodes is not unambiguously known. For the chamber studied here and the EEI impedances for the considered frequency range, we found no significant difference between the lumped and distributed approach (Figure 10), which suggests that the field estimates by the lumped-element approach are reliable. Figure 10 shows also the homogeneity of the electric field.