Predictive Cell Culture Time Evolution Based on Electric Models

Citation: Serrano, J.A.; Pérez, P.; Daza, P.; Huertas, G.; Yúfera, A. Predictive Cell Culture Time Evolution Based on Electric Models. Biosensors 2023,13, 668. https:// doi.org/10.3390/bios13060668 Received: 6 May 2023 Revised: 5 June 2023 Accepted: 6 June 2023 Published: 20 June 2023 Copyright: © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). biosensors Article Predictive Cell Culture Time Evolution Based on Electric Models Juan Alfonso Serrano 1, Pablo Pérez 1,2, Paula Daza 3, Gloria Huertas 1,4 and Alberto Yúfera 1,2,* 1Instituto de Microelectrónica de Sevilla (IMSE-CSIC), Av. Americo Vespuccio 24, 41092 Sevilla, Spain; [email protected] (J.A.S.); [email protected] (P.P.); [email protected] (G.H.) 2 Departamento de Tecnología Electrónica, ETSII, Universidad de Sevilla, Av. Reina Mercedes sn, 41012 Sevilla, Spain 3Departamento de Biología Celular, Facultad de Biología, Universidad de Sevilla, Av. Reina Mercedes sn, 41012 Sevilla, Spain; [email protected] 4Departamento de Electrónica y Electromagnetismo, Facultad de Física, Universidad de Sevilla, Av. Reina Mercedes sn, 41012 Sevilla, Spain *Correspondence: [email protected] Abstract: Obtaining cell concentration measurements from a culture assay by using bioimpedance is a very useful method that can be used to translate impedances to cell concentration values. The purpose of this study was to find a method to obtain the cell concentration values of a given cell culture assay in real time by using an oscillator as the measurement circuit. From a basic cell–electrode model, enhanced models of a cell culture immersed in a saline solution (culture medium) were derived. These models were used as part of a fitting routine to estimate the cell concentration in a cell culture in real time by using the oscillation frequency and amplitude delivered by the measurement circuits proposed by previous authors. Using real experimental data (the frequency and amplitude of oscillations) that were obtained by connecting the cell culture to an oscillator as the load, the fitting routine was simulated, and real-time data of the cell concentration were obtained. These results were compared to concentration data that were obtained by using traditional optical methods for counting. In addition, the error that we obtained was divided and analyzed in two parts: the first part of the experiment (when the few cells were adapting to the culture medium) and the second part of the experiment (when the cells exponentially grew until they completely covered the well). Low error values were obtained during the growth phase of the cell culture (the relevant phase); therefore, the results obtained were considered promising and show that the fitting routine is valid and that the cell concentration can be measured in real time by using an oscillator. Keywords: bioimpedance; cell culture; computer-aided design (CAD); electric model; fractional order (FO); microelectrode; oscillation-based test (OBT) 1. Introduction In recent years, a plethora of scholars have conducted research on monitoring the growth of a cell culture (CC) assay to develop noninvasive, cheap, and robust techniques [ 1 – 9 ]. Biomedical setups have included toxicology assays [ 10 , 11 ], cancer characterization experiments [ 12 – 14 ], biochemical experiments [ 15 ], immune assays [ 16 ], stem cell differentiation protocols [ 17 ], etc., and scholars have sought to use these setups to quantify the number of cells present to characterize a diversity of research objectives and techniques [ 18 – 20 ]. The modeling of a biological sample (BS) allows one to know its electrical behavior and several useful parameters. In many cases, these models are not useful to obtain the required result by themselves. For this reason, the main purposes of this experiment were to enhance the existing models and develop an automatic procedure to obtain the concentration of a CC in real time (RT). The ECIS (electrical cell–substrate impedance sensing) technique is used to sense the electrical response generated on a BS ([ 21 – 23 ]), the CC, when it is excited with an Biosensors 2023,13, 668. https://doi.org/10.3390/bios13060668 https://www.mdpi.com/journal/biosensors Biosensors 2023,13, 668 2 of 15 alternating current or electrical voltage source at several frequencies due to its conductivity properties. The current–voltage relationship, or Ohm’s law, returns an impedance value with real and imaginary parts, the so-called bioimpedance (BI). Because the current applied to the BI must be very low, the ECIS technique requires precise and robust electronics [ 24 ]. Currents below 10 µ A are employed to avoid any damage to the cells, and the applied voltage amplitudes are below 100 mV to correctly bias the electrodes in the linear region [ 5 ]. The CC is immersed in an ionic solution, the culture medium, and the CC cells settle on the substrate of the culture well, where the electrodes are placed. The ideal use of this technique is injecting a signal whose amplitude and frequency are set to their optimum values, i.e., the values that are most sensitive to a change in the BI. From this, the BI value is obtained, but to understand its value, a model of the electrical behavior of an electrode–cell system immersed in an ionic solution must be used. Some scholars have focused on modeling the BI by using solution electrical system equations [ 1 , 5 ] and performing finite elements simulations [25,26] of the whole system, which comprises the CC and electrodes. The authors of reference [ 5 ] modeled the electrochemical behavior of an electrode immersed in an ionic solution. The electrical behavior of an electrode immersed in an ionic solution is modeled as the resistance of the electrode (R ct ) in parallel with a capacitor (C dl ), whose total impedance is Z e (Figure 1a). The spreading resistance (R s ) is the opposition to the current flow in the saline solution that is in contact with the electrode. When the CC grows, it covers the electrode as a cell layer. The cell layer acts as an impedance whose resistive effects are added to the effect of the electrode. By using this model, the real electrical behavior of the cell–electrode (CE) system was reported in [ 27 ]. Figure 1b shows the layout of the electrical model components of a real CC. The electrical model contains the division of the R ct and C dl in two parallel branches, which models the electrode area (A e ) covered by cells (A C ) and not covered by cells (A nC ). In this sense, the electrode is covered by cells (A C ), which represents a measure of the CC time evolution, which is useful to determine the cell number and density. Biosensors 2023, 13, x FOR PEER REVIEW 2 of 16 The ECIS (electrical cell–substrate impedance sensing) technique is used to sense the electrical response generated on a BS ([21–23]), the CC, when it is excited with an alternating current or electrical voltage source at several frequencies due to its conductivity properties. The current–voltage relationship, or Ohm’s law, returns an impedance value with real and imaginary parts, the so-called bioimpedance (BI). Because the current applied to the BI must be very low, the ECIS technique requires precise and robust electronics [24]. Currents below 10 μA are employed to avoid any damage to the cells, and the applied voltage amplitudes are below 100 mV to correctly bias the electrodes in the linear region [5]. The CC is immersed in an ionic solution, the culture medium, and the CC cells settle on the substrate of the culture well, where the electrodes are placed. The ideal use of this technique is injecting a signal whose amplitude and frequency are set to their optimum values, i.e., the values that are most sensitive to a change in the BI. From this, the BI value is obtained, but to understand its value, a model of the electrical behavior of an electrode– cell system immersed in an ionic solution must be used. Some scholars have focused on modeling the BI by using solution electrical system equations [1,5] and performing finite elements simulations [25,26] of the whole system, which comprises the CC and electrodes. The authors of reference [5] modeled the electrochemical behavior of an electrode immersed in an ionic solution. The electrical behavior of an electrode immersed in an ionic solution is modeled as the resistance of the electrode (Rct) in parallel with a capacitor (Cdl), whose total impedance is Ze (Figure 1a). The spreading resistance (Rs) is the opposition to the current flow in the saline solution that is in contact with the electrode. When the CC grows, it covers the electrode as a cell layer. The cell layer acts as an impedance whose resistive effects are added to the effect of the electrode. By using this model, the real electrical behavior of the cell–electrode (CE) system was reported in [27]. Figure 1b shows the layout of the electrical model components of a real CC. The electrical model contains the division of the Rct and Cdl in two parallel branches, which models the electrode area (Ae) covered by cells (AC) and not covered by cells (AnC). In this sense, the electrode is covered by cells (AC), which represents a measure of the CC time evolution, which is useful to determine the cell number and density. Figure 1. (a) Model of electrode immersed in an ionic solution. (b) Basic cell–electrode model. Because real-time cell growth is the main parameter to be measured, the electrical model must include a parameter that acts as a cell-growth indicator. This parameter is the fill factor (ff), which is the percentage of the electrode area covered by cells in %1; i.e., if Ac = Ae, the ff value is 1 (confluence phase), but if Ac = 0, the ff value is 0 (before seeding the cells at t = 0). Thus, Rct1, Cdl1, Rct2, and Cdl2 (the electrode model parameters covered and not covered by the cells, respectively) are defined by the following expressions: 11 22 (1 ) 1 ct ct dl dl ct ct dl dl R R C C ff ff R R C C ff ff = =  − − = =  (1) Figure 1. (a) Model of electrode immersed in an ionic solution. (b) Basic cell–electrode model. Because real-time cell growth is the main parameter to be measured, the electrical model must include a parameter that acts as a cell-growth indicator. This parameter is the fill factor (ff), which is the percentage of the electrode area covered by cells in %1; i.e., if A c =A e , the ff value is 1 (confluence phase), but if A c = 0, the ff value is 0 (before seeding the cells at t = 0). Thus, R ct1 ,C dl1 ,R ct2 , and C dl2 (the electrode model parameters covered and not covered by the cells, respectively) are defined by the following expressions: Rct1=Rct 1−f f Cdl1=Cdl ·(1−f f ) Rct2=Rct f f Cdl2=Cdl ·f f (1) The transfer function (TF), which models the CE impedance, is obtained by solving the circuit in Figure 1b. Equation (2) shows the BI TF, Z CE (s), which is parameterized. It is described by using a parametrized second-order system, where ω0z and Qare its Biosensors 2023,13, 668 3 of 15 natural frequency and quality factor, respectively. The constants k 0 ,k 1 , and k 2 are defined in Equation (2): ZCE(s) = k2·s2+k1·ω0z Q·s+k0·ω2 0z s2+ω0z Q·s+ω2 0z k2=Rsk1=Rs+Rgap·Rct1 2·Rgap+Rct1+Rct2k0=Rs+Rct1·(Rgap+Rct2) Rgap+Rct1+Rct2 ω0z=rRgap+Rct1+Rct2 Rgap·(Cdl ·Rct)2Q=ω0z·Rgap·Cdl·Rct 2·Rgap+Rct1+Rct2 (2) This basic model, the single-electrode well (SEW) model, considers the fact that under the cell, the electrodes are a big, unique electrode. We propose two more models (which are enhancements of the SEW model) of the electrode in Figure 2, which consider the real electrode composed by 10 microelectrodes and 1 large reference electrode: the real-electrode well (REW) and fractional-order well (FO REW) model. Enhanced electrodes with optimum sensitivities can be researched, as shown by the authors of [28]. Biosensors 2023, 13, x FOR PEER REVIEW 3 of 16 The transfer function (TF), which models the CE impedance, is obtained by solving the circuit in Figure 1b. Equation (2) shows the BI TF, ZCE (s), which is parameterized. It is described by using a parametrized second-order system, where ω0z and Q are its natural frequency and quality factor, respectively. The constants k0, k1, and k2 are defined in Equation (2): 22 0 2 1 0 0 22 00 1 1 2 2 1 0 1 2 1 2 12 00 212 () () 2 ( ) 2 zz CE zz gap ct ct gap ct s s s gap ct ct gap ct ct gap ct ct gap dl ct zz gap dl ct gap ct ct k s k s k Q Zs ss Q R R R R R k R k R k R R R R R R R R R R R C R Q R C R R R R     +   +  = +  +   + = = + = +  + + + + + +   = =     + + (2) This basic model, the single-electrode well (SEW) model, considers the fact that under the cell, the electrodes are a big, unique electrode. We propose two more models (which are enhancements of the SEW model) of the electrode in Figure 2, which consider the real electrode composed by 10 microelectrodes and 1 large reference electrode: the real-electrode well (REW) and fractional-order well (FO REW) model. Enhanced electrodes with optimum sensitivities can be researched, as shown by the authors of [28]. Figure 2. 8W10E PET applied biophysics electrodes. Eight separated wells with ten circular biocompatible gold microelectrodes with 250 μm diameter and a large reference electrode with an area larger than each circular microelectrode by approximately 400 squared micrometers. This experiment is based on the experiments and measurements conducted by the authors of previous papers [27]. For these experiments, the measuring circuit was an oscillator. From the oscillation frequency and amplitude, which varies with ff, the cell concentration in real time (RT) can be obtained by using the Barkhausen stability criterion (BSC), which is the mathematical condition that the closed-loop feedback system must fulfill before obtaining sustained oscillations. The objective of the presented experiment is to evaluate the impact of the proposed electric models on the cell cultures and electrodes in the assay and their time evolution in such a way that allows us to predict the cell number and density from the start time of the experiment until its confluence phase. We will then describe the electric models that we employed and the computer algorithms that we developed to create real-time predictions of the cell-culture status. 2. Material and Methods The combination of model fitting, the use of an oscillation-based test (OBT) as a sensor, and the minimization of the cost function (CF) resulted in the ff values and the cell concentration in RT, i.e., during the real growth experiment of a CC assay, being predicted Reference Electrode Sensing Electrodes 8 7 6 5 Reference 1 2 3 4 Figure 2. 8W10E PET applied biophysics electrodes. Eight separated wells with ten circular biocompatible gold microelectrodes with 250 µ m diameter and a large reference electrode with an area larger than each circular microelectrode by approximately 400 squared micrometers. This experiment is based on the experiments and measurements conducted by the authors of previous papers [ 27 ]. For these experiments, the measuring circuit was an oscillator. From the oscillation frequency and amplitude, which varies with ff, the cell concentration in real time (RT) can be obtained by using the Barkhausen stability criterion (BSC), which is the mathematical condition that the closed-loop feedback system must fulfill before obtaining sustained oscillations. The objective of the presented experiment is to evaluate the impact of the proposed electric models on the cell cultures and electrodes in the assay and their time evolution in such a way that allows us to predict the cell number and density from the start time of the experiment until its confluence phase. We will then describe the electric models that we employed and the computer algorithms that we developed to create real-time predictions of the cell-culture status. 2. Material and Methods The combination of model fitting, the use of an oscillation-based test (OBT) as a sensor, and the minimization of the cost function (CF) resulted in the ff values and the cell concentration in RT, i.e., during the real growth experiment of a CC assay, being predicted as the main targets. First, we designed a routine (based on previous work [ 29 ]) in which a real experiment is simulated and the ff value and the cell concentration is predicted at each moment. This routine will be tested by using the REW and FO REW electrical models on the data with three cell lines obtained from the experiment. The first cell line was formed by Chinese hamster ovarian fibroblasts. This cell line was identified as AA8 (purchased from the American type culture collection), and this sample was immersed in McCoy’s medium supplemented with 10% (v/v) fetal calf serum, 2 mM L-glutamine, 50 µg/mL streptomycin, and 50 U/mL penicillin. Biosensors 2023,13, 668 4 of 15 The second and third biological samples that were tested were two mouse neuroblastoma cell lines. The N2a cell line stably expressed a wild-type human amyloid precursor protein, N2a-APP. Both N2A and N2A_APP were provided by Dr. Javier Vitorica from IBiS (Instituto de Biomedicina de Sevilla) Sevilla (Spain). The cells were cultured in a medium consisting of 50% DMEM high glucose (Biowest, Nuaille, France) and 50% Opti-MEM (Gibco, Alcobendas, Spain) supplemented with 10% (v/v) fetal bovine serum (FBS) (Gibco), 2 mM L-glutamine, 50 µ g/mL streptomycin, and 50 U/mL penicillin (Sigma-Aldrich, Madrid, Spain). N2a-APP was also supplemented with 0.4% geneticin (Gibco). All the cell lines were maintained at 37 ◦ C in a humidified atmosphere with 5% CO 2 , and they were routinely subcultured. Different initial numbers of cells (N ini ) were seeded for our experiments: 2500, 5000, and 10,000 cells. The AA8 experiments started with a N ini of 2500 cells in wells 1 and 3, 5000 cells in wells 4 and 5, and 10,000 cells in wells 7 and 8. Moreover, the N2a and N2aAPP experiments began with a N ini of 2500 cells in wells 2 and 6, 5000 cells in wells 3 and 7, and 10,000 cells in wells 4 and 8. The plates contained eight separated wells with ten circular biocompatible gold microelectrodes with a 250 µ m diameter each, which were designed for general-purpose cell culture applications. The electrodes employed in the cell culture assays (8W10E PET), whose electric model was employed for this experiment, were delivered by Applied Biophysics [ 30 ]. These electrodes are fabricated for ECIS equipment and are sold by this company. They are specifically designed for cell cultures and are employed in many assays that are referenced in the bibliography. The electrode size is defined by removing the top isolating mask. The selected size (250 µ m diameter) is related to the common size of cells (1–100 µ m diameter). The sensing area of the electrode is concentrated at the central part of the cultureware, and it obtains the maximum cell sensitivity (to avoid border effects). This layout has been enhanced by the same company, which delivers an IDE configuration (8W10E+ PET) that allows the sensing area of each well to be optimized. An optimized design can be obtained for each cell line, but this is not the purpose of this experiment. The electric models described here can be used for electrical simulation by considering the electrode area as the main parameter to define the final design. The ECIS technique relies on the attachment of the cells to the substrate, which is the electrodes for our experiment. So, a larger electrode area means a higher system sensitivity. This electrode area must be large enough to correctly sample the number of cells. From the datasheet of electrodes, 500–1000 cells can be sensed. This number must be the mean value considering the most common cell-size value, and the cell density must be homogeneous at every cultureware. 2.1. Bioimpedance Modeling 2.1.1. Real-Electrode Well Model (11 Electrodes) The REW model contains the real electrode structure, ten microelectrodes, and one big reference electrode (Figure 2). Figure 3shows a block diagram of the REW model between voltage V CE and the ground. Ten microelectrodes are present, e 1 to e 10 , as well as the reference electrode e r . Each electrode is modeled after the electrode in Figure 1. The only difference is the value of the reference electrode parameters because this depends on the electrode area. To reduce the complexity of the model, the relationship between the microelectrode model parameters and the reference electrode model can be used (R gap is considered to have the same value for microelectrodes and reference electrodes). The relationship between R ct and R ctr (the reference electrode R ct parameter), between C dl and C dlr (the reference electrode C dl parameter), and between R s and R sr (the reference electrode Rsparameter) can be derived from [5] as Rctr Rct =R0 ct Are Ae R0 ct →Rctr =Rct Ae Are [Ω] Cdlr Cdl =C0 dl ·Aer C0 dl ·Ae→Cdlr =Cdl Are Ae[F] Rsr Rs=ρ·√π 4·√Are 4·√Ae ρ·√π→Rsr =RsqAe Are [Ω] (3) Biosensors 2023,13, 668 5 of 15 Biosensors 2023, 13, x FOR PEER REVIEW 5 of 16 sidered to have the same value for microelectrodes and reference electrodes). The relationship between Rct and Rctr (the reference electrode Rct parameter), between Cdl and Cdlr (the reference electrode Cdl parameter), and between Rs and Rsr (the reference electrode Rs parameter) can be derived from [5] as ' ' ' ' [ ] [ ] 4 [ ] 4 ctr ct e e ctr ct ct re ct re dlr dl er re dlr dl dl dl e e e sr e sr s s re re R R A A RR R A R A C C A A C C F C C A A A RA RR RA A   = → =   = → =    = → =   (3) Figure 3. Real-electrode well (REW) model of the well. It has ten microelectrodes (e1 to e10) and one big reference electrode (er). where Are is the reference electrode area and ρ is the electrolyte resistivity of the saline solution. Because the Ae/Are factor is repeated, ke can be defined as ke = Ae/Are. Therefore, the parameters of the electric model of the reference electrode are C dl ctr ct e dlr sr s e e C R R k R R k k =  = =  (4) By using the relationship between the parameters of the microelectrodes and reference electrodes model, Zwell (s) can be defined as the TF of the REW model. Zwell (s) is composed of a gain, one pole, and one zero. The TF is very complex, but it can be simplified when ff→0 (the electrode does not have any cells) and ff→1 (the electrode is full of cells or is in the confluence phase), and one can obtain the pole, zero, and gain expressions for both cases: VCE ZCE(s) e1 ZCE(s) e2 ZCE(s) e10 ZCE(s) er C1 R1Rs C2 R2 Rgap C1r R1r Rsr C2r R2r Rgap Figure 3. Real-electrode well (REW) model of the well. It has ten microelectrodes (e 1 to e 10 ) and one big reference electrode (er). where A re is the reference electrode area and ρ is the electrolyte resistivity of the saline solution. Because the A e /A re factor is repeated, k e can be defined as k e =A e /A re . Therefore, the parameters of the electric model of the reference electrode are Rctr =Rct ·keCdlr =Cdl ke Rsr =Rs·pke(4) By using the relationship between the parameters of the microelectrodes and reference electrodes model, Z well (s) can be defined as the TF of the REW model. Z well (s) is composed of a gain, one pole, and one zero. The TF is very complex, but it can be simplified when ff → 0 (the electrode does not have any cells) and ff → 1 (the electrode is full of cells or is in the confluence phase), and one can obtain the pole, zero, and gain expressions for both cases: pf f →0=−1 Cdl Rct zf f →0=−Rct·(1+10·ke)+Rs·(1+10·√ke) Cdl Rct Rs·(1+10·√ke) pf f →1=−1 Cdl Rct zf f →1=−Rct·(1+10·ke)+Rs·(1+10·√ke)+11·Rgap Cdl Rct·(Rs·(1+10·√ke)+11·Rgap) kf f →0 well =Rct ·1 10 +ke+Rs·1 10 +√ke kf f →0 well =Rct ·1 10 +ke+Rs·1 10 +√ke+11 ·Rgap 10 (5) By using the equations of poles and zeros, the R ct and C dl values when ff → 0 and ff → 1 are derived (because p ff→0 =p ff→1 , only p ff→1 is used in the following expressions to reduce complexity and simplify the model fitting): Rf f →0 ct =−Rs·(1+10·√ke)·(pf f →0−zf f →0) pf f →0·(1+10·ke) Cf f →0 dl =−1+10·ke Rs·(1+10·√ke)·(pf f →0−zf f →0) Rf f →1 ct =−(Rs·(1+10·√ke)+11·Rgap)·(pf f →0−zf f →0) pf f →0·(1+10·ke) Cf f →1 dl =−1+10·ke (Rs·(1+10·√ke)+11·Rgap)·(pf f →0−zf f →0) (6) To increase the fitting, the R s parameter is considered to change with ff, as described in [ 27 ]. Then, R s is split into two terms, R si and ∆ R s , and Rs(k)=Rsi +f f n(k)·∆Rs . The expression means that when ff → 0, R s takes the value of R si , and the ff value on each moment kincreases the Rsvalue. Biosensors 2023,13, 668 6 of 15 2.1.2. Fractional Order Model Fractional order (FO) models are based on the premise that the order of a differential operator can be a noninteger. The differential operator can be defined as the FO constant phase element (CPE), which is used to characterize electrodes for bioimpedance measurements of animal tissue. In combination with reference [ 31 ], the CPE can substitute the C dl term in the REW model to obtain the FO model. The capacitors of the CE model are replaced by the CPE and are described with FO operators: XCdl1=1 sα1Cdl1 XCdl2=1 sα2Cdl2 (7) where α1 and α2 are the FO orders of the reactance XCdl1 and XCdl2 , respectively. Then, as in the previous section, the TF for the cases ff →0 and ff→1 can be obtained: Zf f →0 well (λ1) = kf f →0 well pf f →0 zf f →0 λ1+zf f →0 λ1+pf f →0 Zf f →1 well (λ2) = kf f →1 well pf f →1 zf f →1 λ2+zf f →1 λ2+pf f →1 where λn=sαn0<αn<2 (8) Note that when ff → 0, α2 has no influence on the system behavior, and when ff → 1, α1 has influence on the system behavior. Because changes in the α1 and α2 modify the magnitude slope and the frequency response phase, ∆ R s ([ 27 ]) becomes redundant and can be removed from the FO model. Due to the high complexity of the model (the cross products of α1 and α2 ), the model is implemented in a different way than outlined in the previous section. The REW model is completely implemented to obtain its parameters, whereas the FO models are implemented in a transitional mode. Thus, considering that the pole is constant for any ff value, a transition from zf f →0 to zf f →1 and from kf f →0 well to kf f →1 well is implemented by using ff to change Zf f →0 well (λ1) to Zf f →1 well (λ2) . Then, the implemented model is Zwell(λ1,λ2) = kwell p z λ+z λ+p(9) where p=pf f →0=pf f →1and z= (1−f f )·zf f →0+f f ·zf f →1 kwell = (1−f f )·kf f →0 well +f f ·kf f →1 well λ= (1−f f )·λ1+f f ·λ2 (10) 2.2. Cost Function During the RT estimation, the previously described models utilize the oscillation frequency and amplitude values ( fCE and aCE ) and return the ff and cell concentration values. In the next section, the process of using the model to predict ff in real time is explained, but a parameter (including ff) needs to be obtained from the BSC in some way. The BSC is the mathematical condition that the closed-loop feedback system must fulfill to obtain the sustained oscillations. An oscillator must meet some conditions to obtain self-sustained oscillation, which should include a linear part, G(s=j ω ) in the Laplace domain, and a nonlinear part, N(a osc ,f osc ), where Gis the transfer function of the linear part of the circuit, sis the Laplace operator, jis an imaginary unit, ωosc is the oscillation frequency on rad/s whereby its value ωosc = 2 π f osc , and Nis the linearized model of the nonlinear part of the circuit (an electronic comparator for the present experiment). Note that f osc and a osc are not the same variables as fCE and aCE , but f osc and a osc can be easily estimated by using some internal gains of the oscillator circuit. As the objective of this experiment is not to explain the oscillator circuit in detail, the calculation of the oscillating Biosensors 2023,13, 668 7 of 15 variables will not be explained in detail because its computation is trivial. Then, according to the BSC, the condition that the circuit must meet is f(aosc,fosc)≡1+G(s=jωosc)·N(aosc,fosc) = 0 (11) where f(a osc ,f osc ) is a complex function that depends on f osc and a osc . This function can be rewritten as f(aosc,fosc) = h1(aosc,fosc) + j·h2(aosc,fosc)(12) where h 1 and h 2 are the real and imaginary parts of f(a osc ,f osc ). The main goal of using BSC is obtaining the oscillation parameters. Then, because the oscillation condition function must be equal to 0, in phasorial form, the condition is f(aosc,fosc)≡h∠ϕ=0∠0◦ h=qh2 1+h2 2=0 ϕ=tan−1h2 h1=0◦ (13) where hand ϕ are the module and angle of f(a osc ,f osc ), respectively, and must be equal to 0 and 0 ◦ , respectively, to meet this condition. Equations of the real and imaginary part of the BSC now exist, which must be satisfied to obtain self-sustained oscillations. Using these equations, two parameters of the system could be obtained if all other parameters are known, but in this case, for the ff estimation process, obtaining these parameters is not possible. Therefore, the cost function is minimized to ensure that the BSC is satisfied so that the values of more than two of the model parameters can be obtained for each moment in combination with the fitting routine (which is outlined in the next section). The best way to meet the condition is to use the complex number h(a osc ,f osc ) module. Then, the CF can be defined as h(aosc,fosc)≡qh1(aosc,fosc)2+h2(aosc,fosc)2(14) The fact that h 1 (a osc ,f osc ) and h 2 (a osc ,f osc ) are squared assures that they cannot compensate each other. 2.3. Fitting Routine The key problem with knowing ff at each time during a real experiment is that the oscillation frequency and amplitude ( fCE and aCE ) values are not available when the well is totally covered by cells ( f f →1 ). That is, you can estimate the parameters of the models at the beginning of the experiment when any cells are in the well ( f f →0 ), but certain parameters have no influence on the behavior of the system at this point. These parameters are R gap , zf f →1 , ∆ R s (for the REW model), and α2 (for the FO REW model). The designed routine is as follows: 1. Estimate initial fCE and aCE values : Figure 4shows the block diagram of step 1, where the initial routine is presented in graphic form. During the first hours or days of the experiment, the mean of the last 5 fCE and aCE values is calculated ( fCE and aCE ). As the sampling time (time between measurements) is 1 h, the average of the last 4 h is taken together with the values that were just obtained. After each measurement, after calculating the mean, a check is performed to verify whether the values obtained are greater than the mean of the new measurement plus a margin (k m = 1.005). If this condition is met, as shown in (15), the lowest fCE and aCE are stored as the minimum values. Figure 4also defines the initial R gap value and the value of the constant k m . Note that the index jis the time index and goes from 1 to j max . When calculating fCE and aCE ,jis incremented from 1 until (15) is satisfied. j max is the maximum jvalue, and its value is defined by the number of measurements taken during the real experiment. Biosensors 2023,13, 668 8 of 15 fCE(j)>fCE(j−1)·km&aCE(j)>aCE(j−1)·km(15) Biosensors 2023, 13, x FOR PEER REVIEW 9 of 16 Figure 4. Step 1: block diagram of step 1 of RT simulation, whose target is to find the mean around the 𝑓 𝐶𝐸 and 𝑎𝐶𝐸 minimums. 2. Computation of the initial parameters of the electrical models: Using the minimum 𝑓𝐶𝐸     and 𝑎𝐶𝐸      estimated in the previous step, the initial parameters of the electrical models are fitted. The prediction is performed by using the CF minimization method. For the REW model, the parameters p(ff→1) (whereby p(ff→0) ≈ p(ff→1)), Rsi, and z(ff→0) are calculated, and the Rct(ff→0) and Cdl(ff→0) values can be derived from the parameters by using the two top equations of (6). For the FO REW model, the parameters p(ff→1) (whereby p(ff→0) ≈ p(ff→1)), Rs, z(ff→0), and α1 are calculated, and the Rct(ff→0) and Cdl(ff→0) values can be derived from the parameters by also using (6). The initial parameters that are calculated are the same for all t(j), and therefore, the values are not re-esti- mated during the simulation. The whole process of estimating the initial parameters is illustrated in the Figure 5 block diagram, which starts from the results of step 1 and ends at the beginning of the third and last step. 3. Real-time ff estimation: The last step and the goal of the routine is to predict the parameter ff in RT. Once the initial parameters of the models are obtained (after the previous step), ff is computed for all the previous measurements and the measurements that will be performed until the end of the experiment. Figure 6 describes the whole prediction process. Yes No Simulation starts Step 2 600 1.005 1 gap m R k j = = = 1 ( 1) ( (max(1, 4): ) ( 1) ( (max(1, 4): ) CE CE CE CE jj f j mean f j j a j mean a j j =+ − = − − = − ( ) ( ) ( ) ( 1) & ( ) ( 1) CE CE m CE CE m f j f j k a j a j k  −   −  Figure 4. Step 1 : block diagram of step 1 of RT simulation, whose target is to find the mean around the fCE and aCE minimums. 2. Computation of the initial parameters of the electrical models : Using the minimum fCE and aCE estimated in the previous step, the initial parameters of the electrical models are fitted. The prediction is performed by using the CF minimization method. For the REW model, the parameters p (ff→1) (whereby p (ff→0) ≈ p (ff→1) ), R si , and z (ff→0) are calculated, and the R ct(ff→0) and C dl(ff→0) values can be derived from the parameters by using the two top equations of (6). For the FO REW model, the parameters p (ff→1) (whereby p (ff→0) ≈ p (ff→1) ), R s ,z (ff→0) , and α1 are calculated, and the R ct(ff→0) and C dl(ff→0) values can be derived from the parameters by also using (6). The initial parameters that are calculated are the same for all t(j), and therefore, the values are not re-estimated during the simulation. The whole process of estimating the initial parameters is illustrated in the Figure 5block diagram, which starts from the results of step 1 and ends at the beginning of the third and last step. 3. Real-time ff estimation : The last step and the goal of the routine is to predict the parameter ff in RT. Once the initial parameters of the models are obtained (after the previous step), ff is computed for all the previous measurements and the measurements that will be performed until the end of the experiment. Figure 6describes the whole prediction process. Biosensors 2023,13, 668 9 of 15 Biosensors 2023, 13, x FOR PEER REVIEW 10 of 16 Figure 5. Step 2: block diagram of step 2 of RT simulation, whose target is to find the initial parameters (𝑓𝑓 → 0) by using the 𝑓 𝐶𝐸     and 𝑎𝐶𝐸      values. Figure 6. Step 3: block diagram of step 3 of RT simulation. It describes the simulation process from j = 1 to j = jmax after obtaining the initial parameters. First, the time index j is initialized to start the ff estimation from j = 1 to j = jmax. With the minimization, the model attempts to obtain the values of the parameters for each j measure: ff, Rgap, z(ff→1), and ΔRs (for the REW model) and ff, Rgap, z(ff→1), and α2 (for the FO REW model) by using 𝑓𝐶𝐸(𝑗) and 𝑎𝐶𝐸(𝑗). For this purpose, a loop is used to increment the Step 3 Step 1 REW model REW model or FO REW model CF minimization of REW model ( ) ( ) 00 min , , , CE CE ff si ff h a f p R z →→ → CF minimization of FO REW model ( ) ( ) 0 0 1 min , , , , CE CE ff si ff h a f p R z  →→ → FO REW model Step 2 1 0 mm n =+ = 1nn=+ max nn== max mm== No No REW model REW model or FO REW model CF minimization of REW model FO REW model 0m= CF minimization function 1jj=+ Best fit parameter values 00 22 ( ) (min( ( , ))) ( ) (min( ( , ))) ( ) (min( ( , ))) (REW model) ( ) (min( ( , ))) else ( ) (min( ( , ))) end gap gap ff ff ss ff j ff fval m n R j R fval m n z j z fval m n if R j R fval m n j fval m n  →→ = = =  =  = No max jj== Simulation ends CF minimization function Yes Yes Yes CF minimization of REW model ( ) ( ) 1 min ( ), ( ) ( , ), ( , ), ( , ), ( , ) CE CE gap s ff h a j f j ff m n R m n R m n z m n →   ( ) ( ) 12 min ( ), ( ) ( , ), ( , ), ( , ), ( , ) CE CE gap ff h a j f j ff m n R m n z m n m n  →  0j= Figure 5. Step 2: block diagram of step 2 of RT simulation, whose target is to find the initial parameters ( f f →0)by using the fCE and aCE values. Biosensors 2023, 13, x FOR PEER REVIEW 10 of 16 Figure 5. Step 2: block diagram of step 2 of RT simulation, whose target is to find the initial parameters (𝑓𝑓 → 0) by using the 𝑓 𝐶𝐸     and 𝑎𝐶𝐸      values. Figure 6. Step 3: block diagram of step 3 of RT simulation. It describes the simulation process from j = 1 to j = jmax after obtaining the initial parameters. First, the time index j is initialized to start the ff estimation from j = 1 to j = jmax. With the minimization, the model attempts to obtain the values of the parameters for each j measure: ff, Rgap, z(ff→1), and ΔRs (for the REW model) and ff, Rgap, z(ff→1), and α2 (for the FO REW model) by using 𝑓𝐶𝐸(𝑗) and 𝑎𝐶𝐸(𝑗). For this purpose, a loop is used to increment the Step 3 Step 1 REW model REW model or FO REW model CF minimization of REW model ( ) ( ) 00 min , , , CE CE ff si ff h a f p R z →→ → CF minimization of FO REW model ( ) ( ) 0 0 1 min , , , , CE CE ff si ff h a f p R z  →→ → FO REW model Step 2 1 0 mm n =+ = 1nn=+ max nn== max mm== No No REW model REW model or FO REW model CF minimization of REW model FO REW model 0m= CF minimization function 1jj=+ Best fit parameter values 00 22 ( ) (min( ( , ))) ( ) (min( ( , ))) ( ) (min( ( , ))) (REW model) ( ) (min( ( , ))) else ( ) (min( ( , ))) end gap gap ff ff ss ff j ff fval m n R j R fval m n z j z fval m n if R j R fval m n j fval m n  →→ = = =  =  = No max jj== Simulation ends CF minimization function Yes Yes Yes CF minimization of REW model ( ) ( ) 1 min ( ), ( ) ( , ), ( , ), ( , ), ( , ) CE CE gap s ff h a j f j ff m n R m n R m n z m n →   ( ) ( ) 12 min ( ), ( ) ( , ), ( , ), ( , ), ( , ) CE CE gap ff h a j f j ff m n R m n z m n m n  →  0j= Figure 6. Step 3: block diagram of step 3 of RT simulation. It describes the simulation process from j= 1 to j=jmax after obtaining the initial parameters.