Using a Digital Twin of an Electrical Stimulation Device to Monitor and Control the Electrical Stimulation of Cells in vitro

Electrical stimulation for application in tissue engineering and regenerative medicine has received increasing attention in recent years. A variety of stimulation methods, waveforms and amplitudes have been studied. However, a clear choice of optimal stimulation parameters is still not available and is complicated by ambiguous reporting standards. In order to understand underlying cellular mechanisms affected by the electrical stimulation, the knowledge of the actual prevailing field strength or current density is required. Here, we present a comprehensive digital representation, a digital twin, of a basic electrical stimulation device for the electrical stimulation of cells in vitro. The effect of electrochemical processes at the electrode surface was experimentally characterised and integrated into a numerical model of the electrical stimulation. Uncertainty quantification techniques were used to identify the influence of model uncertainties on relevant observables. Different stimulation protocols were compared and it was assessed if the information contained in the monitored stimulation pulses could be related to the stimulation model. We found that our approach permits to model and simulate the recorded rectangular waveforms such that local electric field strengths become accessible. Moreover, we could predict stimulation voltages and currents reliably. This enabled us to define a controlled stimulation setting and to identify significant temperature changes of the cell culture in the monitored voltage data. Eventually, we give an outlook on how the presented methods can be applied in more complex situations such as the stimulation of hydrogels or tissue in vivo.


STIMULATION CHAMBER
: Images of the stimulation chamber (6-well lid with platinum electrodes) designed for flexible stimulation applications (top view A and bottom view B) as well as a drawing of the device (C).

UNCERTAINTY QUANTIFICATION
A M e a n S D  Figure S3: Comparison of probability distributions sampled from surrogate UQ models for field (A) and current (B). Note that this figure corresponds to the UQ results presented in Figure 3 in the manuscript.

DIRECT CURRENT STIMULATION -CHRONOAMPEROMETRY
We could describe the recorded current response I by a function of the following form: 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 numerical results for the fit parameters are summarised in Table S1.  Table S1. The parameters obtained by fitting the experimental data to Equation (S1). Note that the large value for c, which also has a large uncertainty, for an applied voltage of 1.25 V indicates that this parameters could not be reliably estimated. This could, for example, suggest that the non-faradaic current in this particular case is time-independent because the term exp(−t/c) of Equation (S1) tends to one for c → ∞.
The results for the currents recorded at voltages other than 1 V are shown in Figures S4 and S5.
The faradaic and nonfaradaic contributions to the current are shown separately for each used voltage in Figure S6.  Figure S4: Recorded currents at a fixed voltage of 1.25 V together with results of fit to Equation (S1). The current with reversed polarity is shown in the right panel for comparison.  Figure S5: Recorded currents at a fixed voltage of 1.5 V together with results of fit to Equation (S1). The current with reversed polarity is shown in the right panel for comparison.

1V
1.25V Figure S6: Recorded currents at fixed voltages of 1 V, 1.25 V and 1.5 V together with results of fit to Equation (S1) and the individual contributions of the fitting function. The current with reversed polarity is shown in the right panel for comparison.

FITTING OF IMPEDANCE SPECTRA
We used an equivalent circuit model described in Ragoisha et al. (2010) and shown in Figure S7. If applicable, we used a lead inductance in the equivalent circuit model. The lead inductance was required when we did not directly connect the potentiostat to the chamber. The basic elements of the circuit ( Figure S7) and their characteristic frequency response are shown and explained in Figures S8 and S9. Evidently, the shape of the impedance in the Nyquist plot in Figure 5 in the manuscript cannot be described by a linear superposition of a resistor and a CPE because it cannot be described by a straight line with a fixed slope. Similarly, the impedance computed from the FFT response (shown in Figure 6 of the manuscript) has to contain contributions from an inductance because its phase turns positive at high frequencies.
C dl R ohm R CPE Figure S7: Equivalent circuit model as described in Ragoisha et al. (2010). One part of the interface impedance is modelled by a constant-phase element (CPE) in series with a resistor (R). The other part is modelled by a capacitor (C dl ), which is meant to account for the ionic double layer. The ohmic resistance of the electrolyte (R ohm ) is connected in series with the interface element. To highlight the source of the ohmic resistance, we also refer to it as R medium . If applicable, a lead inductance can be connected in series with the ohmic resistance.  Figure S8: Impedance response of the basic elements of the circuit shown in Figure S7 for a capacitor of 50 mF, a CPE (described by κ(jω) −α ) with κ = 10 and α = 0.5, a resistor of 0.5 Ω and a lead inductance of 100 nH. The values were chosen arbitrarily. This figure is solely intended to highlight the impedance response of the basic elements in the relevant frequency range. The capacitor has a complex-valued impedance with zero real part and negative imaginary part, which appears as vertical line in the Nyquist plot. The CPE has a complex-valued impedance with both non-zero real and imaginary part, which appears as a straight line with a fixed slope starting from the origin in the Nyquist plot. The resistor has a real-valued impedance, which appears as a single point in the Nyquist plot. The inductance has a complex-valued impedance with zero real part and positive imaginary part, which appears as vertical line in the Nyquist plot. The impedance response of a real system is given by linear and parallel combinations of these basic elements.  Figure S7 for a capacitor of 50 mF, a CPE (described by κ(jω) −α ) with κ = 10 and α = 0.5, a resistor of 0.5 Ω and a lead inductance of 100 nH. The values were chosen arbitrarily. This figure is solely intended to highlight the characteristic appearance in Bode plot of the impedance response of the basic elements in the relevant frequency range. The capacitor has a fixed phase of −90 • , the inductance has a fixed phase of 90 • , the resistor has a fixed phase of 0 • and the CPE has a fixed phase of −45 • , which is due to α = 0.5. Aside from the resistor, all elements have frequency-dependent absolute values of the impedance. The impedance response of a real system may contain contributions from these elements, which vary over the measured frequency range. CPEs and capacitors influence mostly the low-frequency range (because the magnitude of their impedance decreases with increasing frequency), while inductances play only a role at higher frequencies (see, for example, Figure 6 of the manuscript). Figure S10: Impedance spectra of KCl solution with 3.5 ml (blue curves), 4 ml (orange curves), and 5 ml (green curves).

Nonlinear EIS
The impedance amplitude and phase at a fixed frequency of 130 Hz for increasing voltage amplitudes are shown in Figure S11.

NUMERICAL COMPUTATION OF STIMULATION WAVEFORMS OF RECTANGULAR TYPE
The complex-valued Fourier series of a signal f is given as where t is time, ω 0 the fundamental angular frequency and c k are the complex-valued Fourier coefficients. The corresponding fundamental frequency is here denoted by f 0 . For a monophasic rectangular pulse of amplitude A and pulse width t p , the Fourier coefficients read Note that in this particular case, the coefficients c k are real-valued. The amplitude of the frequency components A k , which is relevant to estimate the strength of each harmonic, is given by The DC component (c 0 ) is only determined by the amplitude, pulse width and fundamental frequency In practice, the signal can be constructed using because the negative k-values would only affect the imaginary part, which for our application has no relevance.
Using this information, a biphasic rectangular pulse can be easily constructed as a superposition of two monophasic square waves. The square waves have to be time-shifted and their amplitudes need to have opposite signs. A time shift of the signal by t s can be achieved by multiplying Equation (S3) by e −jkω 0 ts . For a biphasic pulse, the time shift has to be equal to the pulse width t p . Writing the signals as a function f (t, A, f 0 , t p , t s ), the signal of the biphasic rectangular pulse reads The amplitudes at the respective frequencies can then be obtained from Equation (S4). Similarly, a biphasic pulse with delay t D can be constructed by Obviously, in both pulses the DC component cancels out because it only depends on amplitude and pulse width.
The presented approach is numerically very simple to implement and permits more flexibility than deriving the analytical presentation for each waveform. A drawback of the Fourier series approach are large oscillations that occur at the jumps of the square waves (the so-called Gibbs phenomenon). We attenuated the oscillations using the Lanczos sigma factor in the summation (Weisstein, 2021). Furthermore, we used 2500 harmonics (i.e., including harmonics up to a frequency of 325 kHz) to compute the signals unless otherwise stated. This covers the relevant part of the frequency spectrum for all considered signals (Figures S12-S14). Using this number of harmonics reduced oscillations caused by high-frequency contributions. This was only evident for current-controlled stimulation, where high-frequency oscillations occurred (Figures S15). They were most likely caused by the increased lead inductance. Such features were not measured. Reasons for that could be a low-pass filter of the stimulator or of the oscilloscope, which was used for recording. Importantly, the general signal shape of the predicted and measured signal agreed well.   Figure S15: Voltage response for a biphasic pulse in current-controlled regime (6.5 mA amplitude) with a pulse width of 200 µs (A) and 600 µs (B) for a single electrode pair. The experimental data is compared to the prediction based on the impedance model (computed using 5000 harmonics). For comparison, the expected voltage is shown: it was computed using the impedance measured by EIS and not inferred from the measured pulses. The overshooting in the prediction data arises due to the Fourier series approach and the high-frequency contribution of the lead inductance. It can be removed by using less harmonics (see Figure 8 in the manuscript, where 2500 harmonics were used while not changing the actual waveform).    Figure S20: Comparison of stimulation voltages (peak-to-peak) for current-controlled stimulation (6.5 mA amplitude, 130 Hz, 60 µs) when the system was used directly after handling without thermal equilibration (blue) or after thermal equilibration took place (i.e., at 37 • C) (orange). The stimulation was performed with three (out of six) filled wells in series (without thermal equilibration) or six filled wells (with thermal equilibration). For the sake of comparison, the results for the thermally equilibrated system were divided by two to account for the additional wells. Because we could show that the wells are equal from an electrochemical point of view, this approach is feasible and does not introduce an error. When the thermally unequilibrated system was used, three 6-well plates, containing fresh cell culture medium, were stimulated successively either for 30 minutes, 60 minutes, or 120 minutes including short handling breaks in between. During these handling breaks, the lid was cleaned and kept at room temperature and the ambient temperature in the incubator dropped by a few degrees Celsius. The stimulation of the thermally unequilibrated system was completed after 4 hours. The stimulation of the thermally equilibrated system was continuously applied for 12 hours. Voltage and current responses (Figures 7, 8, S15, S16, S17, S18, S19)

Entities
Activities prior characterisation  Figure S21: Provenance graph showing entites and activities. We show three major activities as well as the entity "Model". The relations (i.e., arrows) indicate which entities are being generated or used by which activities.