Modeling and simulation of non-linear bioelectronic systems applied to cell culture assays

Modeling and simulation of non-linear bioelectronic systems applied to cell culture assays Juan Alfonso Serrano Viseas Supervised by: Dr. Alberto Yúfera García Dra. Gloria Huertas Sánchez December 2021 2 I would like to dedicate this thesis to my parents, my brother and all my family, for their support throughout this journey. To Alicia, for encouraging me at every step of the way. To all my friends, for each and every one of their phrases in the best and worst situations. To my colleagues at the Seville Institute of Microelectronics (IMSE-CNM), without whose contribution this thesis would not be possible. To my tutors, for their guidance during these years. 3 Acknowledgements This thesis would not have been possible without the help of many people. Each and every one of them has provided both technical and support in other ways. If I neglect to mention another name, just remember that I have not forgotten your support and that you have done your bit to make this thesis possible. This work has been possible thanks to funding from the University of Seville and the Spanish government with the Mixcell project (REF: TEC2013-46242-C3-1-P) and Istent project (REF: RTI2018-093512-B-C21). The first persons to be mentioned are Prof. Alberto Yúfera García and Prof. Gloria Huertas Sánchez. They have guided me during these years, providing me with encouragement, technical support and their experience for the successful completion of this thesis. Thank you both for the excellent support you have given me, I am very proud and it has been a pleasure to have you as my tutors. Thanks to my colleagues at Seville Institute of Microelectronics (IMSE-CNM), especially to the PA34 members for helping me in some technical details. I will never forget his friendship and support to follow and complete this journey. I would also like to thank Prof. Antonio José Ginés Arteaga. During almost all this period, and especially during one year, he has shared and transmitted much of his knowledge and experience. I am very grateful for your help and I am very proud to have worked with you. Prof. Paula Daza and Prof. María Martín of the Cell Biology Department in the Seville University, for their essential contribution to obtain the required data without which this work could not have been completed. Also, to all the staff of the Institute of Biomedicine of Seville (IBiS), where extra data were obtained from cell cultures and tissues. Thanks to all of you for helping to obtain the biological data that have made this work possible. My parents, Alfonso and Concepción, and my aunt, Encarna, have always encouraged me and helped me to choose the right decisions at this stage and in all the previous stages of my life. Also, to my brother, Antonio, who helped me with his ideas in some technical aspects to continue in some moments of this thesis. Thank you very much, for everything. Finally, thank you Alicia for caring for me in every bad moment, for helping me smile when my ideas didn't work, and for being my daily support during almost all this stage of my life. 4 5 Abstract Monitoring the properties of biological samples (BS) is expensive in terms of time consumption and cost in resources and human effort. There are monitoring methods using different techniques, depending on the type of biological sample (cells, tissues, blood, etc.). In this thesis, we improve the accuracy of the ECIS (electrical cell-substrate impedance spectroscopy) technique of bioimpedance (BI) measurement of a monolayer cell-culture (CC), which measures the electrical response of a CC when an alternating current is applied to it at several frequencies. Specifically, this dissertation is focused on the modeling of the cell-electrode (CE) block, and the real-time monitoring and acquisition of the cell concentration in a CC assay experiment. In addition, using as a database the CE model obtained in the modeling stage, a toolbox has been built to perform efficient electrical real-time simulations with Ngspice, launching these simulations from Matlab. Herein, the BI measurement is applied indirectly. Instead of injecting a signal, the CE block is connected to an electronic oscillator, which fulfills the Barkhausen Stability Criterion (BSC) to ensure that self-maintained and self-sustained oscillations are generated. The technique is known as Oscillation Based Test (OBT) [1]. Instead of measuring the changes between the output and input signal, the oscillation frequency and amplitude are acquired and, using the BSC, the parameters of the CE model used can be obtained on the fly. As can be seen, this technique is much more powerful than injecting a signal, since instead of obtaining the BI for a given frequency, the whole CE electrical model of the CE block is obtained. The CE model used is based on previous work, but some improvements are introduced to increase its accuracy. The data base for this work are real measurements made by the research group in CC assays with three different cell lines using an OBT circuit. From these data, the parameters of the CE model block are successfully obtained. The modeling technique is tested on some variations of the CE model, reaching better results by increasing its complexity, making the model closer to reality or introducing Fractional Oden (FO) elements. Data from real experiments, and the best variants of the electrical model are used to build a simulator of a CC assay experiment. The simulator calculates model parameters and cell concentration in real-time (without taking into account future measurements) using the minimization of a cost function (CF). The minimization of an appropriate CF ensures that the oscillation requirements are satisfied and that the obtained CE model parameter values are consistent with the theoretical values. The acquired results are very satisfactory since the simulation of a real-time experiment demonstrates that the technique of minimizing a CF can be used to obtain the cell concentration in real-time. As a result, cell concentration data are attained whose trend and values present a relatively low error when compared to the cell concentration achived by traditional optical cell counting methods. Electrical simulations of an electronic circuit are very useful during the design and testing process prior to the build of the circuit. The model parameters obtained during 6 the simulation of a CC assay experiment, in addition to showing the feasibility of the technique, can be used in electrical simulations of the OBT circuit. For future improvements of the OBT measurement circuit, which are discussed at the end of this dissertation, electrical simulations must be performed with realistic data to ensure that the measurement circuit will work robustly. To perform such simulations, a toolbox has been built in Matlab, which performs electrical simulations using Ngspice (open-source Spice simulator) in an efficient way. This toolbox is applied to the oscillator circuit simulation successfully, being able to perform multiple simulations, varying the CE model automatically, in an efficient way. In addition, it can be applied to any electronic circuit to launch electrical simulations from Matlab. 7 Table of contents List of figures………………………………………………………………………9 List of tables...…………………………………………………………………….13 Glossary of abbreviations……………...……………………………………….17 Chapter 1 Introduction ...................................................................................... 19 1.1 Biological samples ................................................................................................................ 20 1.2 Electrical modeling and applications ................................................................................ 22 1.3 Challenges of cell-culture assay modeling ....................................................................... 24 1.4 Data Acquisition solutions .................................................................................................. 26 1.5 Objectives ............................................................................................................................... 27 Chapter 2 Electrical modeling of biological samples .................................. 31 2.1 Cell-culture assay modeling ................................................................................................ 32 2.1.1 Cell-culture assay growth ............................................................................................... 33 2.1.2 Cell-electrode model ........................................................................................................ 34 2.1.3 Measurement circuit theory ............................................................................................ 39 2.1.3.1 Oscillation-based test ............................................................................................... 40 2.1.3.1.1 Describing function method ............................................................................ 40 2.1.3.1.2 Barkhausen stability criterion .......................................................................... 43 2.1.4 OBT measurement circuit ............................................................................................... 44 2.1.4.1 Band-pass filter ......................................................................................................... 44 2.1.4.2 Comparator ................................................................................................................ 46 2.1.4.2.1 Pre-amplifier stage ............................................................................................ 46 2.1.4.2.2 Comparator stage .............................................................................................. 48 2.1.4.3 Current source ........................................................................................................... 49 2.1.5 Empirical cell-electrode model ....................................................................................... 50 2.1.5.1 Cell lines specifications ............................................................................................ 50 2.1.5.2 Single-electrode well model .................................................................................... 51 2.1.5.2.1 Bode fitting method .......................................................................................... 52 2.1.5.2.2 Empirical fitting method using oscillation parameters ................................ 55 2.1.5.3 Real-electrode well model (11 electrodes) ............................................................. 63 2.1.5.4 Fractional order model ............................................................................................. 69 2.1.5.5 Model fitting of different cellular lines .................................................................. 74 2.1.5.5.1 Fitting results of N2a and N2aAPP cell lines ................................................. 74 8 2.1.5.5.2 Results comparison between cellular lines .................................................... 85 2.2 Conclusions ............................................................................................................................ 89 Chapter 3 Simulations of real experiments of biological systems ........... 91 3.1 Real-time estimation of cell-concentration growth ........................................................ 91 3.1.1 Fitting routine ................................................................................................................... 92 3.1.2 Real-time simulation results ........................................................................................... 96 3.1.2.1 AA8 results ................................................................................................................ 96 3.1.2.2 N2a results ............................................................................................................... 102 3.1.2.3 N2aAPP cell line results ......................................................................................... 107 3.1.2.4 Comparison of cell line results .............................................................................. 111 3.2 Electrical simulation of the Bio-Oscillator ..................................................................... 112 3.2.1 Electrical simulator and its integration in Matlab ..................................................... 113 3.2.2 Electrical simulation using REW electric model ........................................................ 114 3.2.2.1 AA8 ........................................................................................................................... 116 3.2.2.2 N2a ............................................................................................................................ 117 3.2.2.3 N2aAPP .................................................................................................................... 117 3.3 Conclusions .......................................................................................................................... 118 Chapter 4 Conclusions and future works .................................................... 121 References…...…………………………………………………………………...127 Appendix A: Scientific publications…………………………………………133 Appendix B: Ngspice toolbox code……………………...……………………135 9 List of figures Fig. 1.1 Growth curve of a CC. There are six different zones which are labeled in the figure. Variable x is the cell concentration on each moment. ....................................................... 21 Fig. 1.2 Dielectric constant (ε) and the conductivity (σ) change depends on frequency. There are three different dispersion zones related to the frequency ( α, β and γ). Figure from [10]. .......................................................................................................................................... 22 Fig. 1.3 Model of an electrode immersed in an saline solution, reported by [17]. Cdl is the interfacial capacitance, Rct is the charge transfer resistance, Rs represent the solution (or spreading resistance) and Cw and Rw are the Warburg elements which are related to diffusion effects. ................................................................................................................ 24 Fig. 1.4 Model of an electrode immersed in a saline solution covered by cells. Z(ω) is the impedance of the electrode immersed in a saline solution, ω is the angular frequency, Ae is the electrode area, and Acell is the area covered by cells. .......................................... 25 Fig. 1.5 Basic BI measurement circuit. Current flow through the BS, voltage between the edges of the BS is measured, and the Ohm`s Law is used to calculate the BI or Z(ω) . .................................................................................................................................................. 26 Fig. 1.6 Cell concentration estimation project. A CC is connected as a load to an OBT circuit. The changes in the oscillation parameters (fosc and aosc) are used by an algorithm to estimate the cell concentration of the CC at each moment. ............................................. 28 Fig. 2.1 Each CC property can be used to know other ones, since some mathematical relation could exist. .............................................................................................................................. 32 Fig. 2.2 8W10E PET Applied BioPhysics electrodes. Eight separated wells with ten circular biocompatible gold microelectrodes of 250 μm diameter, and a large reference electrode with an area approximately 400 larger than each circular microelectrode... 34 Fig. 2.3 Electrical model of an electrode immersed in an ionic solution impedance without the presence of cells. .................................................................................................................... 35 Fig. 2.4 Double-layer interphase diagram, where double layer effect is divided in two terms: electrode covered by cells and electrode non-covered by cells. ...................................... 37 Fig. 2.5 Double-layer interphase model. a) CE diagram seen from the top, which shows the electrode area covered by cells (ZC), the electrode area not covered by cells (ZnC), the reference electrode and the current flow lines flowing through electrode and cells to reference electrode. This is a diagram for illustration purpose, electrodes have not the real size seen at Fig. 2.2. b) Equivalent electrical circuit that models the electrical behavior of the CE system. ................................................................................................... 38 Fig. 2.6 Basic block diagram of a self-sustained oscillator. ............................................................ 40 Fig. 2.7 Non-linear linear feedback system where there are two different parts, the linear (𝐺(𝑗 · 𝜔𝑜𝑠𝑐)) and the non-linear part ( 𝑁(𝑎𝑜𝑠𝑐, 𝑓𝑜𝑠𝑐)). ..................................................... 41 Fig. 2.8 Comparator hysteresis (a), comparator output (b) and comparator input (c). ............. 42 16 17 Glossary of abbreviations AC Alternative Current BSC Barkhausen Stability Criterion BI Bioimpedance BPF Band-Pass Filter BD Bode Diagram BS Biological Sample CE Cell-Electrode CC Cell-culture CPE Constant Phase Element CF Cost Function DF Descriptive Function DE Differential Equation DC Direct Current ECIS Electrical Cell-substrate Impedance Spectroscopy EC Environmental Conditions ES Electrical Simulation FE Finite Elements FO Fractional Order HPF High-Pass Filter IO Integer Order LPF Low-Pass Filter NA Network Analyzer OBT Oscillation-Based Test OCF Oscillation Condition Function OHP Outer Helmholtz Plane PI Proportional-Integral REW Real-Electrode Well RT Real-Time SUT Sample Under Test SEW Single-Electrode Well TF Transfer Function 18 19 Chapter 1 Introduction Throughout history, humanity has tried to understand everything that is around it. It started with relatively simple-to-understand elements until today, where it is possible to treat complex diseases, design and fabricate electronic circuits, etc. The study of nature, and all the microorganisms that reside in it, is a subject that, although great progress has been achieved, is still in the process of research. Biologists from all over the world are currently involved in the study of cells, tissues, blood studies, etc. Many tools used in this kind of research have progressed a lot, such as electron microscopes, incubators, and other tools that are used today. However, there are still processes such as optical cell counting, cell differentiation analysis, disease detection, etc. that are still too "rudimentary" compared to the enormous advances seen in other research areas. In the area of cell biology, specifically in Cell Culture (CC) assays studies, CC monitoring methods require a very large amount of time and resources. All tasks and experiments performed with CC require one or more people to monitor the status of the culture every time. This means, for example, that the researcher must count the number of cells 1 by 1 in each culture (the traditional method of counting is explained in detail later), check their morphology, size, etc. This is done, for evident reasons, using a microscope. In order to observe a cell under a microscope, it must be dead and have been stained with a colorant. Once the cells are killed and stained, this culture can no longer be monitored since, as mentioned above, the cells (the CC) have died. Researchers, in order to know the state of a CC assay, must perform a complex and expensive process. To take a single measurement of the cell concentration present in a well, they must take a CC assay, kill and stain the cells, and then perform an optical count through a microscope. The optical (or traditional) counting with the microscope is statistical, several samples (or zones) of the culture are counted and extrapolated to the whole CC using the mean of the taken measurements. This process is long and, moreover, is expensive, since for each measurement that biologist want to perform, they must "waste" a CC assay. Considering the present level of technology, it is coherent to think that it is possible to develop an easier and more efficient method to perform all these tasks. In recent years, some methods to solve this problem have emerged, such as image processing to calculate the number of cells in a CC assay [2], [3], flow cytometry (a sophisticated instrument measuring multiple physical characteristics of a single cell) [4], or bioimpedance (BI) measurement [5], [6]. The BI is the resistance to the current flow produced by any biological sample (BS). There are many studies on BI, its modeling and applications, but there are also many unsolved problems, such as the translation of BI into useful parameters for the field of biology, such as cell concentration, differentiation level, composition of a cell, etc. This work focuses on the use of BI to try to solve or go a step further in solving this problem these points. Concretely, the main goal is to find a method to measure the cell concentration of a CC assay using the BI. 20 The structure of the present document begins with Chapter 1, which contains an introduction, where the current state of the art of the agents with which it has been working (CC assays, BI modeling methods, measurement circuits, data acquisition and processing, and simulation) will be explained. The current state of the art of BS and its electrical modeling is reviewed. After that, the challenges of BS modeling is presented. The fourth point is to choose and compare the different ways of data acquisition when taking electrical measurements of biological samples. This is a very complex field, because there are many types of sensors and acquisition methods. Finally, the starting point and objectives of this thesis are presented. In Chapter 2, the model of the chosen measurement circuit and the possible models that can be used to predict the behavior of a CC assay are studied. Details are given about how the chosen measurement circuit, an oscillator, works. On the other hand, different ways of processing the data (bode diagrams, Barkhausen criterion, etc.) are studied. Finally, the parameters of three variants of the same model are fitted using data from a CC growth experiment, and the results are compared with each other and with data obtained in the traditional experimental way (optical counting with microscope). Then, in Chapter 3, two types of simulations applicable to BI measurement are presented and tested. The first one is the simulation of a cell growth experiment in which the cell concentration is estimated in real-time. To perform real-time estimation electric models of the CC must be used, which is the reason why these models have been studied in Chapter 2. The second one is the electrical simulation (ES) and its use for the design of BI measurement circuits, exposing an efficient method to perform electrical simulations applied to this work. In the future, these electrical simulations will be very useful when an improvement of the measurement circuit is required, since the real-time simulation data is used to predict the behavior of the circuit. Finally, Chapter 4 presents the conclusions reached in this study and possible future work that might be derived from it. 1.1 Biological samples A BS is a sample of biological material from humans, animals or plants, such as blood, cells, tissue, DNA, etc. Each type of sample has its own composition, morphology and size, which can be measured in different ways, which are defined as properties of a BS. The present work focuses on the study of CC assays, which possess properties of their own that can be quantified. For example, cell differentiation is a process by which cells are transformed into more complex cell types or even tissues. The differentiation level is a fundamental property of a CC assay when performing a cell differentiation study [7] and there are many techniques for sensing and acting on tissues [8]. Another property, and one of the targets of this thesis, is the growth of a CC assay. The growth or proliferation of a cell culture is the process by which the number of cells and / or the size of cells in a CC assay grows. This process can be of the cell division type, in which a cell divides into two smaller ones; cell growth, which occurs when a cell increases in 21 size; growth and cell division, in which the cell increases its size to then divide into two cells of the same size as the original cell. There are many works devoted to the study of cell proliferation, such as [9], which details the process of cell counting by optical density (a less accurate technique than direct counting). Fig. 1.1 shows the growth curve of a CC, which has been reproduced from [9], which also details the six different zones (labeled with numbers in the figure) of the CC growth. 1. Lag phase: The cells adapt to the new culture medium and the new environmental conditions. 2. Acceleration phase: Some structural components of the cells increase sequentially and the cell concentration in the well begins to increase. 3. Exponential phase: The growth rate reaches its maximum and the cell concentration increases quickly. 4. Deceleration phase: the cells face fewer growth factors that enhance growth, so the growth rate of the concentration decreases. 5. Steady-state phase: The birth rate is equal to the mortality rate, so the growth rate of the CC concentration is negligible. 6. Apoptosis (death phase): The mortality rate exceeds the birth rate, so the growth rate is negative. This means that the cell concentration decreases. As previously explained, cell counting by traditional methods (or direct counting) is very laborious and expensive in terms of time and cost. The researchers must spend time preparing several CCs in order to count, for example, each day, the number of cells in 1 2 3 4 5 6 Time ln (x) Fig. 1.1 Growth curve of a CC. There are six different zones which are labeled in the figure. Variable x is the cell concentration on each moment. 22 that CC. If a fairly accurate growth curve is required (Fig. 1.1), it is necessary to kill, stain and count cells from several similar CC cell cultures (seeded at the same time) on the same day. Even then there may be errors with respect to what would be a growth curve measured each day on the same CC assay, since the counting is based on obtaining measurements from several CC and obtaining the statistical mean. As can be deduced, this process is not optimal, since it requires the expenditure of many materials (cells, petri dishes, etc...), and a lot of time performing the whole process. 1.2 Electrical modeling and applications As mentioned, this work focuses, among other topics, on the measurement of BI, a property that all BS exhibit. The BI property is based on the electrical properties reported in [10] in 1957 by Schwan, the dielectric constant (ε) and the conductivity (σ). These two properties are frequency dependent as shown in Fig. 1.2 (obtained from [10]) and three frequency dependent dispersion zones are observed: α, β and γ. Each BS, depending on its morphology and composition, has a different BI. The BI can be translated into one of the properties of the BS, thus making an indirect measurement of the desired property. For example, in [11], the measurement of BI, specifically conductivity, is used to determine the blood glucose level. The study [12] describes a method to estimate interstitial fluid volume and peripheral blood volume in hemodialysis patients using whole body and calf bioimpedance spectroscopy techniques. On the other hand, for CC assays there are a lot of researches that study the Fig. 1.2 Dielectric constant (ε) and the conductivity (σ) change depends on frequency. There are three different dispersion zones related to the frequency ( α, β and γ). Figure from [10]. 23 relationship of CC BI and some of its properties. In [13] the BI of some CCs is tracked to find relationships between BI and cell differentiation. There are also many studies focused on research for a robust, reliable and accurate non-invasive technique to study and estimate cell growth in a CC [5], [14]–[17]. In addition, cell number quantification has also been targeted for various research purposes and techniques, such as toxicity assays [18], cancer characterization experiments [19], [20], biochemical [21], immuneassays [22], stem cell differentiation protocols [23], etc. All of these studies have at least one thing in common, they all seek to use the BI to indirectly quantify one or more properties of a BS. To translate the BI into the desired property, electrical models are used. Electrical models predict the behavior of a given material or object in the presence of electric current, i.e., they indicate how the material or object will react against the flow of current. Furthermore, depending on the complexity level of the model, they will be composed of a certain number of parameters. These parameters are usually linked to some property of the material or object, such as its density, size, shape, etc. In the case of BS, the same applies. When the electric current flows through a BS, it faces an impedance, which depends on the type of BS (blood, cells, tissue, etc.) and its properties and characteristics (size, shape, composition, concentration of some chemical substance, etc.). These parameters are the key to obtain and quantify, indirectly, the properties of a BS. The Bioimpedance-based measurement approach, also called ECIS (electrical cellsubstrate impedance spectroscopy) measures the response of a CC when an alternating current is applied to it at several frequencies. Impedance spectroscopy or ECIS technique is used because it is a reliable and robust non-invasive technique to estimate and study cell growth on cell-culture assays. This technique requires robust electrical models of the cell-electrode (CE) block and accurate electronic circuits to extract the signals to be measured [24]. Note that this technique is based on the idea that in a CC assay, cells are seeded in a well filled with culture medium and adhere to the bottom of the well, thus forming the CE block through which the electric current flows. BI monitoring and modeling has been previously studied in several works [5], [17], [25], [26], [27]. In [5] and [27] the authors solve the electrical equations that model the electrical behavior of the BS under study. Research of [25], [26] use Finite-Element simulations of the whole system (cell-electrode-culture medium system) to extract the relationships between the conductivity, morphology and electrical behavior of the system. Many of the studies mentioned above are the base for this paper, in particular the works [17], [25], [28]. In [17] the model of an electrode immersed in a saline solution (without the presence of cells or other BS) is studied and derived. In [28] a motility experiment is performed with mammalian cells and determines the impedance of a cellcovered electrode and a non-cell-covered electrode. The work [25] performs finite element simulation to find the value of some elements of the electrical models used by the previously described works and, in addition, introduces a new parameter. 24 1.3 Challenges of cell-culture assay modeling The CE model used in this work is based on the model of an electrode immersed in a saline solution (culture medium) [17], including in the model the presence of cells [19]. About this model, promising test have been reported, due to the good results obtained in the cited works. In [17] it is explained how to calculate some of the parameters of the model (Fig. 1.3), which is shown in 2.1.2. The diffusion-related elements (Rw and Cw, named Warburg resistance) are calculated theoretically by [29]. These elements, for the frequencies used in this study, are removed because their effect on the model is negligible. The charge transfer resistance (Rct) depends on exchange current density (J0) and the thermal voltage (Vt), which are parameters that can vary easily, because their theoretical determination depends on the properties of the electrode and the saline solution, as well as temperature and other variables. In the computation of Rct there can be errors due to the number of parameters and their variability, and the same occurs in the calculation of the interfacial capacitance (Cdl), because it depends on much more parameters with similar or greater variability. Finally, as for the spreading resistance (Rs), its estimation is relatively more accurate and easier than the previous parameters, since it depends on the electrode geometry and its resistivity. However, when cells are involved, Rs is more complex to estimate and can even change, as reported in [25] and [30]. The model of an electrode immersed in a saline solution is very complex and presents a high variability. The calculation of the parameters (2.1.2), although relatively accurate, is complex, since there are many variables that can deviate the value of any of them. One of the main challenges of the present work is to find experimental values, or ranges of values, of these Cdl Rs Cw Rw Rct Metal Electrolyte Fig. 1.3 Model of an electrode immersed in an saline solution, reported by [17]. Cdl is the interfacial capacitance, Rct is the charge transfer resistance, Rs represent the solution (or spreading resistance) and Cw and Rw are the Warburg elements which are related to diffusion effects. 25 parameters. The theoretically estimated parameters are used as a starting point to know the approximate behavior of the model and design the experiments accordingly. An effect not included in [17], is the gap resistance (Rgap), which is found by [25] using Finite-Element simulations. The model, included cells and its effects, like Rgap, is shown in Fig. 1.4, where these effects are shown. Z(ω) is the impedance of the electrode immersed in a saline solution, ω is the angular frequency, Ae is the area of the electrode and Acell is the area of the cells. Z(ω) is divided in two parts, the covered by cells and the uncovered one. In addition, the Rgap effect is added after the impedance of the cell covered electrode. Rgap represents the effect of the current flowing between the cells attached to the bottom of the well and the well itself, and is highly dependent on the space (gap) between cell and well, which can takes values of 15-150 nm as explored in [25], [31]. In [25] Rgap values of 105 Ω are considered, while in [31], the resistance assumes values of 90 kΩ. The model of Fig. 1.4 contains the theoretical behavior of an electrode covered by cells and immersed in a saline solution when an electric current is flowing through it. This model is composed of different parameters that model each of the parts of the system, except for one. Cells are considered to be an infinite resistance, i.e., no electric current flows through them. In other works, the cell impedance is considered as a capacitor [32], or as a capacitor and a resistor [33]. In the present work it has not been considered since the impedance division model of the electrode immersed in a saline solution, combined with Rgap and Rs, seems to predict well the behavior of the system. The parameters of the model are difficult to estimate theoretically, due to the fact that the size and composition of the elements that are part of the system are complex to accurately calculate. In addition, the model contains a considerable number of variables, which must be computed and fitted to the data of a real experiment. Finally, biological elements usually exhibit nonlinear behavior, i.e., they cannot be fully modeled using conventional passive elements. There are techniques, such as the fractional-order (FO) models, which can be applied to the BS modeling. In [34], FO models are used to Rs Rgap Z(ω)/Acell Z(ω)/(Ae-Acell) Fig. 1.4 Model of an electrode immersed in a saline solution covered by cells. Z(ω) is the impedance of the electrode immersed in a saline solution, ω is the angular frequency, Ae is the electrode area, and Acell is the area covered by cells. 32 • Differential equation (DE): A DE is a mathematical equation that relates a function to its derivatives. If one focused on the solution of a time-dependent DE, the relation between u and y can be gotten. A DE could be linear or nonlinear. If it is linear, it can be solved by Laplace transform, i.e., transforming it in a TF. If it is non-linear, there is some solving methods like linearizing in a working point, numerical solving, etc. In this work an electronic oscillator is proposed to monitoring the growth of a CC assay. The CC assay is connected as a load to the oscillator, which produces changes in oscillation frequency and amplitude. The oscillation frequency and amplitude are used, in combination with electrical models of the CC, to predict the growth of the CC assay. This chapter describes the oscillator circuit and the electrical model of the CC assay on which this work is based, as well as all the modifications of the model, which will be tested throughout the chapter. Finally, the model parameters will be fitted to real experimental data performed with three different cell lines and the results achieved will be discussed. 2.1 Cell-culture assay modeling In recent years, a lot of research has been focused on growth monitoring of a cellculture assay for the development of a non-invasive, cheap and robust technique [5], [14]–[17]. Biomedical setups such as toxicology assays [18], cancer characterization experiments [19], [20], biochemical [21], immune-assays [22], stem cells differentiation protocols [23], etc., seek to quantify the number of cells for characterizing a diversity of research objectives and techniques. As was previously introduced, the modeling of a BS allows to know its electrical behavior and several useful parameters. Cell-culture Properties • Cell number • Bioimpedance • Duplication rate • Electrical Model • Size • Shape • . Fig. 2.1 Each CC property can be used to know other ones, since some mathematical relation could exist. 33 ECIS (electrical cell-substrate impedance sensing) technique senses the electrical response generated on a BS ([28], [49]), the CC, when it is excited with an alternating current (AC) electrical source, voltage, or current, at several frequencies, as consequence of its conductivity properties. Since the current applied to the BI must be very low, the ECIS technique requires precise and robust electronics [24]. Usually, the CC is immersed in an ionic solution, the culture medium, and the cells of the CC settle on the substrate of the culture well, where the electrodes are placed. The ideal use of this technique would be to inject a signal whose amplitude and frequency are set at their optimum values, i.e., the values most sensitive to a change in BS. The current-voltage relationship, or Ohm's law, returns an impedance value with a real part and an imaginary part, the so-called BI, which is defined as: BI BI BI Z R j X= +  (2.3) In this way, the value of the BI is obtained, but to understand its value, a model of the electrical behavior of an electrode immersed in an ionic solution must be used. There are some works focused on modeling BI by solve electrical system equations [5], [17] or performing finite elements (FE) simulations [25], [26] of the whole system, which comprises CC and electrodes. Then, the main goal is to obtain a cell-electrode (CE) model which throws up information about the current state of the CC, i.e., about the CC growth. 2.1.1 Cell-culture assay growth CC assay tests are the best way to known the behavior of cells in specific situations. These tests give data about the response of cells in any isolated scenario, i.e., CC assay tests allow to isolate the effect of a certain situation in cells without the effects of another situation. These "situations" could be cell differentiation scenario, toxicology test, unusual cell behavior characterization, etc. For example, in a cell differentiation test, differentiation data could be obtained without the "perturbation" of another effects in cells. In that way, the cell growth is a basic parameter to characterize a CC. CC tests must be performing in a controlled environment, i.e., the environmental conditions (EC) must be properly set for cell proliferation. Tests are held in a CC reactor, or "incubator", in which the researcher can chose the EC. Concretely, the data used in this work are obtained from CC assay experiment whose EC are a temperature of 37 °C and 5 % CO2 in a very humid atmosphere (90 %). CC assay is placed in a Petri dish or in another well type. Cells take some hours (number of hours varies depending on cell line) to adhere to the substrate in the bottom of the Petri dish [9]. Each certain time period, cells are detached from the bottom and rise to perform cell division (mitosis phase), i.e., cell divides into two identical cells and identical to the cell that originated them. In that way, when there are enough cells, the CC builds a layer over the bottom of the Petri dish. Cell division is repeated until all the bottom is covered by cells, obtaining a monolayer. 34 2.1.2 Cell-electrode model ECIS technique is based on the Double-layer interphase explored in [5], [17], [25], [47], among others, and also in the behavior of a CC assay. An electrode immersed in an ionic solution generates an electric-field on the electrode surface. For this purpose, Applied Biophysics [46] commercial electrodes are used. The employed electrodes (model 8W10E PET) contain eight separated wells with ten circular biocompatible gold microelectrodes of 250 μm diameter (Fig. 2.2), which are labeled in Fig. 2.2. Each well also contains a reference electrode (reference electrodes of each well are connected to each other) which has an area approximately 400 larger than each circular electrode. Reference [17] models the electrochemical behavior of an electrode immersed in an ionic solution. This idea is displayed in Fig. 2.3; the electrical model of the impedance of an electrode immersed in an ionic solution (CC medium). That is, the electrical behavior of an electrode immersed in an ionic solution is modeled as a resistance (Rct) in parallel with a capacitor (Cdl), whose total impedance is Ze. In [17], [47] the physical meaning of each component of the electrical model (Fig. 2.3) is explained in detail, so in this work each parameter will be explained and its theoretical estimation: • Charge transfer resistance (Rct) Rct, depends on exchange current density (J0) and the thermal voltage (Vt). J0 is the absolute value of current per area in the electrode, and it depends on the electrode properties and ionic solution. According to [36], referenced by [17] and [47], for a gold micro-electrode immersed in a saline solution, the value is J0 = 3.98·10-6 A/cm2. Thermal voltage is estimated by: 26.71 t kT V mV q  == (2.4) Reference Electrode Sensing Electrodes 8 7 6 5 Reference 1 2 3 4 Fig. 2.2 8W10E PET Applied BioPhysics electrodes. Eight separated wells with ten circular biocompatible gold microelectrodes of 250 μm diameter, and a large reference electrode with an area approximately 400 larger than each circular microelectrode. 35 where T is the temperature in kelvin (K) set at 310.15 K (37 °C), q is the electron charge and it is set to 1.60218·10-19 C, and k is the Boltzmann’s constant that has a value of 1.38·10-23 J/K. The charge transfer resistance per area is calculated by (2.5). Then, the expression to estimate Rct of a gold microelectrode of area Ae immersed in an ionic solution is (2.6). '2 0 671.21 t ct V R G m J  = =  (2.5) ' [ ] ct ct e R RA = (2.6) • Double layer capacitance (Cdl) Over the electrode immersed in an ionic solution, a charged capacitive layer is created. This layer was modeled by Stern [50], basis in works of Helmholtz [51], and Gouy [52] and Chapman [53]. Helmholtz capacitance (CH) in terms of capacity per unit area is defined by: 2 0 [ / ] r H OHP C F m d   = (2.7) where 𝜀0 is the vacuum permittivity (8.85419·10-12 F/m), εr is the relative permittivity of the solution (value of 6.8), and dOHP is the distance between electrode surface and Outer Helmholtz Plane (OHP) ([17], [54]) whose value is 8.25·10-10 m. For these values, CH is 0.073 F/m2. The second capacitance included in that model is the Gouy capacitor (CG), and can be estimated[17] by: 2 0cosh( ) [ / ] 2 ro G Dt zV C F m LV   =  (2.8) Cdl Rct Rs Ze Fig. 2.3 Electrical model of an electrode immersed in an ionic solution impedance without the presence of cells. 36 where z is the valence of the ions present in the reaction of creation of the double layer (z = 1), Vo is the maximum electrode potential (provided by [46]) whose value is 100 mV [47], and LD is the Debye Length, models the drop in potential along its length [17]. LD is defined by: 002 [ ] 2 rt D V Lm n z q   =   (2.9) n0 is the concentration of ions in the electrolyte whose value, for a saline solution, is n0 = 140 mol/m3 [47]. q is the electron charge, q = 1.60218·10-19 C. Then, using this data LD has a value of 2.44·10-10 m and CG = 0.821 F/m2. The total capacitance per unit area (𝐶𝑑𝑙 ′) is the series of Helmholtz capacitance (CH) and Gouy-Chapman capacitance (CG). Cdl of a gold microelectrode of area Ae immersed in an ionic solution can be calculated by (2.11). '2 ' 1 1 1 0.067 / dl dl H G C F m C C C = + → = (2.10) ' [ ] dl dl e C C A F= (2.11) • Spreading resistance (Rs) Spreading resistance is the opposition to current flow through saline solution in contact with electrode. This resistance is generally determined by (2.12), but (2.13) is the equation to calculate Rs for a circular electrode [17]. [ ] s L RA   = (2.12) = [ ] 44 s ee RrA     =  (2.13) where ρ is the electrolyte resistivity, whose value for a saline solution of NaCl is 5·105 Ω·μm, L is the length, A is the cross-sectional area (cm2) of the solution through which the current passes and re is the electrode radius (re = 125 μm). The theoretical value of Rs for a circular electrode of radius re is 1 kΩ. • Gap resistance (Rgap) The last parameter of the basis electrical CE model is Rgap, which is shown in Fig. 2.5. This resistor models the resistive gap region through which the current flows laterally, and it was modeled in [25] using finite elements simulations. Rgap depends on adhesion distance of the cell to the electrode, which varies between 15 nm and 150 nm [31]. In this reference, a value for Rgap of 75 kΩ is considered, but in [30], [55] values near to 1 kΩ was obtained. In the present work, Rgap is estimated from real experimental data When CC growths, it covers the electrode as a cell layer. The cell layer, act as an impedance whose effects is added to the effect of the electrode. In that way, the real 37 electrical behavior of the CE system has been reported in [30], [55]–[58]. Fig. 2.4 and Fig. 2.5 shows the model of CE covered by cells. First shows the layout of the electrical model components in a real CC, while in the second figure illustrates the schematic of the electrical model as seen from above the CC. The electrical model presents the division of Rct and Cdl in two parallel branches, modeling the electrode area covered by cells (AC) and not covered (Ae - Ac) by cells, equation (2.14). That division is shaped by two impedances. () ee C nC C e C ZZ ZZ A A A == − (2.14) where ZC is the impedance of the circuit covered by cells and ZnC is the impedance of the circuit non-covered by cells. Since cell growth is the main parameter to measure, the electrical model must include a parameter that informs as an indicator of cell growth. This parameter is the fill-factor (ff), which is the percentage of electrode area covered by cells in %1, i.e., if AC = Ae the value of ff is 1, but if AC = 0 the value of ff is 0. In that way R1, C1, R2 and C2 are defined by next equations: 11 (1 ) 1 ct ct dl dl R R C C ff ff = =  − − (2.15) 22 ct ct dl dl R R C C ff ff = =  (2.16) where Rct1 and Cdl1 are the contribution of Rct and Cdl in the area not covered by the cells, and Rct2 and Cdl2 are the contribution of Rct and Cdl in the area covered by the cells. The TF which models the CE impedance is obtained by solving the circuit of Fig. 2.5. Equation (2.17) shows the BI TF, which has two poles and two zeros. Poles and zeros are directly dependent on ff (among other things), which means that the number of cells covering the electrode influences the BI model. Equation (2.17) has been parameterized because the equation is too big. Electrode Double Layer Cell Cdl1 Rct1 Cdl2 Rct2 Rgap RsReference Electrode Culture medium Fig. 2.4 Double-layer interphase diagram, where double layer effect is divided in two terms: electrode covered by cells and electrode non-covered by cells. 38 22 0 2 1 0 0 22 00 () zz CE zz k s k s k Q Zs ss Q    +   +  = +  + (2.17) where k2, k1, k0, ωn and Q are the model parameters. These parameters are defined below as: 2s kR= (2.18) 1 1 12 2 gap ct s gap ct ct RR kR R R R  =+  + + (2.19) 12 0 12 () ct gap ct s gap ct ct R R R kRR R R + =+ ++ (2.20) 12 02 () gap ct ct z gap dl ct R R R R C R  ++ = (2.21) 0 12 2 gap dl ct z gap ct ct R C R QR R R   =  + + (2.22) Cdl1 Rct1 Rs Cdl2 Rct2 Rgap ZC ZEC ZnC b) Reference electrode Electrode a) ZnC ZC ICE ICE V1V2 ΔV = V1-V2 Fig. 2.5 Double-layer interphase model. a) CE diagram seen from the top, which shows the electrode area covered by cells (ZC), the electrode area not covered by cells (ZnC), the reference electrode and the current flow lines flowing through electrode and cells to reference electrode. This is a diagram for illustration purpose, electrodes have not the real size seen at Fig. 2.2. b) Equivalent electrical circuit that models the electrical behavior of the CE system. 39 The described model (Fig. 2.5) is the main model proposed in this work, but in next sections will be discussed and improved, for the sake of completeness. 2.1.3 Measurement circuit theory Traditionally, the BI measurement is performed by applying AC current signals of different frequencies [31]. Depending on the frequency and the cell type, each impedance measure is different. By applying a signal, a complex measured impedance is obtained. The impedance value, without a model to translate it to another property, can only be used during the experiment to see the trend of a certain property which is known to vary when impedance changes. In that way, when experiment is over, if the relationship between impedance and the property to be measured is unknown, the typical approximation between the two properties is to find the relationship between the property to be measured and the relative impedance, for example, based on the full range of the observed impedance along the experiment or its initial value. Equations (2.23) show these examples. Where i is the time index of the experiment and its range is 1 ≤ 𝑖 ≤ 𝑘 (where k is the i value where impedance is maximum), 𝑟𝑍(𝑖) is the relative impedance in the i moment, Z(i) is the impedance in i moment, Z(1) is the first impedance measurement, and Z(k) is the maximum and last impedance measurement. These equations show the typical relative impedances used to obtain information from a BS under test: ( ) (1) ( ) (1) ( ) ( ) (1) ( ) (1) ZZ Z i Z Z i Z r i r i Z Z k Z −− == − (2.23) Due to ECIS technique is used to perform an estimation of change in one or more properties of a BS, and this change is related to an impedance change, it is hard to find a relationship between an absolute property value of a CC assay (for example, the cell growth in a CC assay) and a relative property value of a CC assay, the impedance relative variation. By known the impedance in a certain frequency, the model parameters cannot be obtained (the only available parameter is the impedance at this frequency). If the model parameters are not estimated, a lot of information could be ignored or lost. To obtain the value of each model parameter, the measured data must be more than the impedance at a certain frequency. The Oscillation-based test is a good option for this purpose. This measurement technique applied to ECIS is based on an indirect measurement of the BI collecting data of the oscillation parameters, i. e., amplitude (aosc) and frequency (fosc). If the BI is connected to an oscillator, as a load, aosc and fosc depends on the BI features. If BI changes during an experiment the oscillation parameters will change as well. 40 2.1.3.1 Oscillation-based test A non-linear system could produce oscillations of constant amplitude and frequency under some conditions, which are called self-excited oscillations or limit cycles. Van der Pol oscillator proposed one of the first equations to analyze these systems. A self–starting and self-sustained oscillator is a circuit which has two main parts (Fig. 2.6), linear part and non-linear part, and it has to meet some requirements: 1. There is only one non-linear element. 2. Non-linear element is time invariant. 3. Non-linear element is symmetrical. 4. Linear part must have a low-pass filter behavior. If these conditions are met by the circuit, there would exist oscillations at a constant amplitude (aosc) and frequency (fosc). Oscillation parameters (aosc and fosc) depend on linear and non-linear parts of the entire system. These parameters could be estimated by the Barkhausen stability criterion (BSC) [1] (it will be explained later), but to apply it, the non-linear element must be linearized. The Descriptive Function (DF) method is employed to obtain the linearized response of a periodic non-linear component from the truncated Fourier Transform. The model of the system added to these two tools is the way to know the behavior of a non-linear oscillator. 2.1.3.1.1 Describing function method The DF is a method based on the frequency response of the system which ignores the higher harmonics. DF method replaces the non-linear system by the first component of the Fourier series. As mentioned above, if the linear part of the system behaves as a low- y(t) u(t) x(t) ν(t) Non-linear Linear Fig. 2.6 Basic block diagram of a self-sustained oscillator. 41 pass filter or a band pass-filter, the decision to truncate the Fourier series is a reasonable approach. First, it is considered that the input of non-linear element is a sine of amplitude aosc and angular frequency of ω. By making the variable change 𝜃 = 𝜔 · 𝑡, where t is the time variable, the input of non-linear element is 𝑦(𝑎𝑜𝑠𝑐, 𝜃)= 𝑎𝑜𝑠𝑐 · sin 𝜃. Fig. 2.7 shows the non-linear feedback system block diagram, where m is defined by: ( ) ( ( , )) osc m t f y a  = (2.24) The non-linear element DF, N, can be presented by using the Fourier coefficients, in complex notation, as: 12 ( ) ( ) ( ) osc osc osc N a f a j f a= +  (2.25) The next steps and considerations to estimate the DF are based on section 2.2 of [59]. Since the non-linear component is symmetric, the offset Fourier component is equal to 0. f1(aosc) and f2(aosc) coefficients are the normalized Fourier coefficients. 11 12 ( ) ( ) osc osc osc osc ab f a f a aa == (2.26) where a1 and b1 are the Fourier coefficients. Then the normalized coefficients could be estimated by 2 10 1 ( ) [ sin( )]sin( ) osc osc osc f a f a d a     = (2.27) 2 20 1 ( ) [ sin( )]cos( ) osc osc osc f a f a d a     = (2.28) The DF estimation can be applied to the non-linear component of the OBT circuit of this work, an inverter comparator with hysteresis. Fig. 2.8 presents the comparator hysteresis in a), the comparator output in b) and the comparator input in c). These figures show that there are four different zones in terms of θ. The comparator output takes the value M or -M depending on θ, where M is the comparator output according to the input. These regions are defined in the expression (2.30). The angle θ1 has a value of: 11 1 () sin osc y a   − =  (2.29) N(aosc , fosc) y(aosc ,ϑ)G(j·ωosc) m(t) + - x(t)r(t) Fig. 2.7 Non-linear linear feedback system where there are two different parts, the linear (𝐺(𝑗 · 𝜔𝑜𝑠𝑐)) and the non-linear part ( 𝑁(𝑎𝑜𝑠𝑐, 𝑓 𝑜𝑠𝑐)). 48 As BPF filters of pre-amplifier stage circuits are implemented using OPA4227UA opamps. TF of the two filters in series is obtained to analyze the entire system in Matlab® [60] to compare the accuracy of the model with the Spice filters [61]. In Fig. 2.14 are presented the BDs of the filter union. The mathematical model is very accurate since it has a low error compared to Spice electrical simulation. 2.1.4.2.2 Comparator stage The comparator stage includes the comparator, the hysteresis passive elements and the diodes to limit the output amplitude. Fig. 2.15 shows the comparator stage. Comparator is implemented in the circuit using an integrated circuit LMC7211 [62]. The comparator has a response time of 4 μs (at 5 V of supply voltage). Ro is the output resistor. Resistors RH and RF, and capacitor CF compose the hysteresis circuit, which set the thresholds to 20 mV . The last part of the comparator stage is the diode net, where D1 and D2 limit the negative output voltage and D3 and D4 limit the positive output voltage. 101102103104105 10 20 30 Magnitude [V/V] 101102103104105 Frequency [Hz] -500 -400 -300 -200 Phase [º] Matlab Spice Fig. 2.14 High-pass filter and low-pass filter BD. Blue line represents Matlab® simulation, and the red line represents the electrical simulation in Ngspice. + - RH CF RF Ro Vi Vo D1 D2 D3 D4 Fig. 2.15 Comparator stage, which includes hysteresis and voltage limit diodes. 49 The Spice transient simulation of the circuit in Fig. 2.15 can be compared with the expected signals (Fig. 2.8). For this purpose, in Fig. 2.16, the output versus input of the comparator is shown. The comparator has a response time of 6 s  (near of expected in datasheet from LMC7211), and its thresholds voltages are near to ±20 mV (observed in simulations). There is a small difference between the positive and negative threshold, so the response of the comparator is somewhat asymmetrical. Considering the work frequencies of the oscillator (10 Hz to 100 kHz), this is not a significant limit. 2.1.4.3 Current source Current source is responsible for ensuring that the current supplied to the CE system is adequate. If the current is too large it could burn the cells and electrodes and, on the other hand, if it is too low, the measurements could be wrong as well. For that reason, a current source independent of the load has been used. The current will be not dependent on the impedance changes due to cell growth. Fig. 2.17 shows the current source input. Taking into account the relationships R3=R1/2 and R4=R2/2, by analyzing the circuit, the current expression in CE is obtained: 2 15 CE i R IV RR − = (2.44) where ICE is the current that flows through BI. From this expression, the voltage gain of the circuit (gain between Vi and Vo) can be obtained: 22 1 5 1 5 ( ) ( ) ( ) CE CE o i CS R Z s R Z s V V H s R R R R −− = → = Resistor values are chosen to obtain a gain current of ICE = -5·10-6Vi, which means that the current in CE system is in the order of ICE = 7 μA (safe value for cells). -0.1 -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 0.1 -1.5 -1 -0.5 0 0.5 1 1.5 Vo [V] Vi [V] Fig. 2.16 Comparator output presented against input (comparator negative pin). 50 2.1.5 Empirical cell-electrode model In this section the CE model is obtained the from empirical data. 𝑍𝐶𝐸(𝑠) (2.17) parameters are estimated using experimental data obtained from the test performed with some cell lines for the sake of illustration [55], [63]. In next subsections, cell lines specifications and the method to obtain the empirical models are shown. 2.1.5.1 Cell lines specifications For this work, three cell lines have been used. The origin of cells, the culture medium, the growth environment conditions and the number of cells seeded is explained below. The first biological Sample Under Test (SUT) is formed by Chinese hamster ovarian fibroblasts: This cell line is identified as AA8 (American Type Culture Collection). This sample is immersed in McCoy’s medium supplemented with 10 % (v/v) foetal calf serum; 2 mM L-glutamine, 50 μg/ml streptomycin and 50 U/ml penicillin. The second and third biological samples under test are two Mouse neuroblastoma cell lines. The N2a cell line and its cell line stably expressing wildtype human amyloid precursor protein, N2a-APP. Cells were cultured in medium consisting of 50% DMEM High glucose (Biowest) and 50% Opti-MEM (Gibco) supplemented with 10% (v/v) foetal bovine serum (FBS) (Gibco), 2 mM L-glutamine, 50 μg/ml streptomycin and 50 U/ml penicillin (Sigma-Aldrich). N2a-APP was also supplemented with 0.4% Geneticin (Gibco). All cell lines were maintained at 37° C in a humidified atmosphere with 5% CO2 and they were routinely subcultured. Different initial number of cells (Nini) was seeded for our experiments: 2500, 5000 and 10000 cells. The AA8 experiments start with a Nini of Op-Amp + - R2 R1 R3 R4 Op-Amp + - R5 VCE ZCE(s) Bioimpedance Vi ICE Fig. 2.17 Current source with CE connected as load. Current in CE system is independent of BI. 51 2500 cells in wells 1 and 3, 5000 cells in wells 4 and 5, and 10000 cells in wells 7 and 8. On the other hand, the N2a and N2aAPP experiments start with a Nini of 2500 cells in wells 2 and 6, 5000 cells in wells 3 and 7, and 10000 cells in wells 4 and 8. Numbers of the wells were defined in Fig. 2.2. 2.1.5.2 Single-electrode well model Single-electrode well (SEW) model is referred to estimate the empirical model by considering that the cells are on top of a single large electrode, and the effect of the reference electrode (grounded) is negligible. The first step to understand the changes in the model when the CC assay growth is to know to position of the poles and zeros fo the system in the important points of the experiments. These ones are when the experiment is started (𝑓𝑓 → 0) and when the well is full (𝑓𝑓 → 1). Poles (2.45), zeros (2.46) and gain (2.47) of the biological system can be obtained from (2.17), and are: 12 (1 ) 1 (1 ) gap ct dl ct gap dl ct R ff ff R pp C R R ff ff C R −+ − == − (2.45) 2 1 1,2 2 2 2 2 1 22 2 (1 ) 4 2 2 ct s ct gap gap s gap s dl ct gap s ct gap gap s gap s s R R n R R ff R R ff R R ff zC R R R ff ff n R R ff R R ff R R ff R  + + − =− = + − + (2.46) 1 11 CE s ct ct gap kR ff RR Rff =+ − − + (2.47) As can be observed the expressions are very large, but if ff is set to 0 or 1, poles, zeros and gain are greatly reduced: 0 00 () 1 ff ct s ff ff CE ct s dl ct dl ct s RR p z k R R C R C R R → →→ −+ − = = = + (2.48) 1 11 () 1 () ct s gap ff ff ff CE ct s gap dl ct dl ct s gap R R R p z k R R R C R C R R R → →→ − + + − = = = + + + (2.49) Poles and zeros for 𝑓𝑓 → 0 and 𝑓𝑓 → 1 provide key information since the model changes from a two-pole two-zero model to a one-pole one-zero model when the well is empty and full. The limit models, for 𝑓𝑓 → 0 and 𝑓𝑓 → 1, are defined by (2.48) and (2.49) respectively. The limit models are obtained applying the limit when ff tends to 0 and when tends to 1: 52 00 00 00 () ff ff ff ff CE CE ff ff p s z Z s k z s p →→ →→ →→ + =+ (2.50) 11 11 11 () ff ff ff ff CE CE ff ff p s z Z s k z s p →→ →→ →→ + =+ (2.51) In the next sections, two different fitting methods using the experimental BD and the oscillation parameters will be detailed. 2.1.5.2.1 Bode fitting method The first proposed fitting method is to use an experimental bode analysis. For that method, AA8 cells in the amounts indicated in section 2.1.5.1. Cells were seeded in 8W10E PET Applied BioPhysics electrodes. Two wells per initial cell concentration were seeded, i.e., 2x5000 cells (wells 1 and 3) and 2x10000 cells (wells 7 and 8). CC was growing for 6 days. Each day at the same time a frequency analysis was performed using a network analyzer (NA), the HP 3589A Spectrum Network Analyzer , to estimate the BD of each well. The measurement circuit is shown in Fig. 2.18. Vi is connected to the NA output and Vo is connected to the NA input. Fig. 2.19 displays the BD of the wells 7 and 8 (10000 cells at the beginning of the experiment) the first day of the experiment and the day when ff is near to 1 (day 5). From the BDs the cut frequency of the pole and the zero can be obtained: 00 11 Well 7: 24.06 Hz 18.844 Well 8: 22.73 Hz 49.124 ff ff ff ff p z kHz p z kHz →→ →→ == == (2.52) Rin Op-Amp + - Vi ZCE(s) Bioimpedance Vo Network Analyzer Output Input Fig. 2.18 Measurement circuit used to obtain the real BD of the CE system in the network analyzer. 53 Due to poles when 𝑓𝑓 → 0 and 𝑓𝑓 → 1 have similar expressions ((2.48) and (2.49)) and its values are very close to each other, the value of 𝑝𝑓𝑓→0 is taken to the estimation. At this point there are 3 equations (𝑝𝑓𝑓→0, 𝑧𝑓𝑓→0 and 𝑧𝑓𝑓→1) and 4 variables (Rct, Cdl, Rs and Rgap). But the limit of the TFs (2.50) and (2.51) when the frequency (𝑠 = 𝑗 · 𝜔) tends to 0 and ∞ can be used to estimate some parameters and provide the following equations: 0 0 00 0 0 00 ( 0) () lim lim ff ff CE ct s sff ff ff CE s sff sz Z s R R sp sz Z s R sp → → →→ → → →→ + → = = + + + →  = = + (2.53) a) b) ff = 0 ff = 1 Fig. 2.19 BDs of wells 7 and 8. a) Data of the first day of the experiment ( 𝑓𝑓 → 0). b) Data of the fifth day of experiment ( 𝑓𝑓 → 1) 54 1 1 01 1 1 01 ( 0) () lim lim ff ff CE ct s gap sff ff ff CE s gap sff sz Z s R R R sp sz Z s R R sp → → →→ → → →→ + → = = + + + + →  = = + + (2.54) Note that 𝑍𝐶𝐸 𝑓𝑓→0(𝑠 → 0) is the gain 𝑘𝐶𝐸 𝑓𝑓→0 and 𝑍𝐶𝐸 𝑓𝑓→1(𝑠 → 0) is the gain 𝑘𝐶𝐸 𝑓𝑓→1. If the measurement circuit is taken into account (Fig. 2.18), and the magnitude values at the limit frequencies of the BDs are obtained, the resulting expressions are: ( ) ( ) 27.8 020 27.05 020 ( 0) 10 ( ) 10 dB ff ct s CE in dB ff s CE in RR Zs R R Zs R → − → + → = = →  = = (2.55) ( ) ( ) 24 120 26.8 120 ( 0) 10 ( ) 10 dB ct s gap ff CE in dB s gap ff CE in R R R Zs R RR Zs R → − → ++ → = = + →  = = (2.56) Using equations (2.48) and 𝑍𝐶𝐸 𝑓𝑓→0(𝑠 → ∞) in (2.55), model parameters can be obtained for 𝑓𝑓 → 0: 0 0 0 382.14 17.31 488.53 ff ff ff ct dl s R k C nF R → → → =  = =  (2.57) Note that it is not required to estimate 𝑅𝑔𝑎𝑝 𝑓𝑓→0 because it does not affect 𝑍𝐶𝐸 𝑓𝑓→0(𝑠). This is an important point, Rgap has no influence when 𝑓𝑓 → 0, and in following sections it will be seen that it has a low effect for low ff values. Fig. 2.20 shows the comparison between the BD of empirical model with 𝑓𝑓 → 0 and the BD of the real CE system. The biggest errors appear for low and high frequencies, in the extremes. 101102103104105106 -20 0 20 Magnitude [dB] Empirical bode model Experimental data 101102103104105106 Frequency [Hz] 100 150 200 Phase [º] Fig. 2.20 BD of empirical model (blue) when 𝑓𝑓 → 0. It is compared to the BD of the real CE system (red). 55 On the other hand, using equations (2.49) and 𝑍𝐶𝐸 𝑓𝑓→1(𝑠 → ∞) in (2.56), models parameters can be obtained for 𝑓𝑓 → 1: 11 11 393.29 16.81 488.53 R 14.26 ff ff ct dl ff ff s gap R k C nF R →→ →→ =  = =  =  (2.58) 𝑅𝑠 𝑓𝑓→0 has also been used for 𝑓𝑓 → 1 because there is no way to obtain values of 𝑅𝑠 𝑓𝑓→1, 𝑅𝑐𝑡 𝑓𝑓→1 and 𝑅𝑔𝑎𝑝 𝑓𝑓→1 greater than 0. Then, 𝑍𝐶𝐸 𝑓𝑓→1(𝑠 → ∞) from (2.56) has been used to estimate 𝑅𝑔𝑎𝑝 𝑓𝑓→1. The value of 𝑅𝑔𝑎𝑝 𝑓𝑓→1 seems wrong, since [25], [55] reported higher Rgap values. Fig. 2.21 shows the BD of the model when 𝑓𝑓 → 1. As expected, model cannot establish the gain and phase value because in that point of the experiment, the CE system is highly non-linear. Another type of model or fitting method can improve the accuracy of this model. 2.1.5.2.2 Empirical fitting method using oscillation parameters OBT circuit oscillates at a frequency and an amplitude which depend on the load 𝑍𝐶𝐸(𝑠) (Fig. 2.9). The amplitude used to calculate the DF (aosc) in the circuit is the amplitude of the signal VLo in Fig. 2.13, but the amplitude that will be measured to characterize the CE system is the amplitude of the signal VCE from Fig. 2.17, from now called 𝑎𝐶𝐸. Both amplitudes, 𝑎𝐶𝐸 and aosc, are related each other by the gains of the preamplifier stage of the comparator (2.1.4.2.1). Therefore, in this section, oscillation amplitude and frequency (oscillation parameters) are employed to estimate the 𝑍𝐶𝐸(𝑠) parameters. Oscillation parameters are obtained from real experiments where CE system is connected to OBT circuit as a load (Fig. 2.9). These experiments are performed by seeded an initial number of cells (Nini) in wells of Fig. 2.2 (different Nini on each well). CC assay grows for a week, from seeding to 101102103104105106 -20 0 20 Magnitude [dB] Empirical bode model Experimental data 101102103104105106 Frequency [Hz] 100 150 200 Phase [º] Fig. 2.21 BD of empirical model (blue) when 𝑓𝑓 → 1. It is compared to the BD of the real CE system (red). 56 confluence, in which the circuit takes a measurement of the oscillation parameters every hour. Oscillation frequency (𝑓𝐶𝐸) and amplitude (𝑎𝐶𝐸) are used in an algorithm to estimate the 𝑍𝐶𝐸(𝑠) parameters. BSC can be applied to OBT circuit to obtain the oscillation parameters. In this case, there is the inverse problem, oscillation parameters are known and can be used to calculate 𝑍𝐶𝐸(𝑠) model parameters. Equations (2.48) and (2.49) describe the position of the poles and zeros when 𝑓𝑓 → 0 and 𝑓𝑓 → 1. Since Rct >> Rs and Rct >> Rgap, the zeros equations can be approximated as: 01 11 () ff ff dl s dl s gap zz C R C R R →→ −− == + (2.59) BSC applied to the OBT system gets the following equation (where Laplace variable is 𝑠 = 𝑗𝜔𝑜𝑠𝑐): ( , ) 1 ( ) ( ) ( ) ( ) ( , ) osc osc BP K HPF LPF CS osc osc h a f H s H H s H s H s N a f +      (2.60) The expression in equation (2.60) must be equal to zero in both module and phase. Then, two equations can be derived from BSC: 1 2 ( , ) ( ( , )) 0 ( , ) ( ( , )) 0 osc osc osc osc osc osc osc osc h a f real h a f h a f imag h a f == == (2.61) where real and imag are the real and imaginary part of ℎ(𝑎𝑜𝑠𝑐, 𝑓 𝑜𝑠𝑐). The algorithm used to estimate oscillation parameters is simple, and it is divided in two steps: 1. Initial estimation. Frequency and amplitude values and equations (2.61) are used to find the initial parameters. The main goal of this step is to obtain 𝑅𝑐𝑡 𝑓𝑓→0, 𝐶𝑑𝑙 𝑓𝑓→0, 𝑅𝑠 𝑓𝑓→0, 𝑝𝑓𝑓→0 and 𝑧𝑓𝑓→0 parameters. Rgap has no influence in system behavior at this point (𝑓𝑓 → 0), so it is defined by 600 Ω. 𝑝𝑓𝑓→0 is calculated because his value is required for the second step. To minimize the number of parameters to be estimated, equation (2.48) and the approach of 𝑧𝑓𝑓→0 in (2.59) is used to define 𝑅𝑐𝑡 𝑓𝑓→0and 𝐶𝑑𝑙 𝑓𝑓→0 dependent on 𝑅𝑠 𝑓𝑓→0, 𝑝𝑓𝑓→0 and 𝑧𝑓𝑓→0. Then, 𝑅𝑐𝑡 𝑓𝑓→0and 𝐶𝑑𝑙 𝑓𝑓→0are defined by (2.62) and (2.63), and used in 𝑍𝐶𝐸(𝑠) in (2.60). 0 0 0 1 2 ff ct ff ff dl RpC  → → → = (2.62) 0 0 0 1 2 ff dl ff ff s CzR  → → → = (2.63) 2. Final estimation. When wells are fully covered by cells, the oscillation parameters can be used to estimate the CE model parameters when 𝑓𝑓 → 1. In the same way that step 1, 𝑅𝑐𝑡 𝑓𝑓→1and 𝐶𝑑𝑙 𝑓𝑓→1are defined by (2.64) and (2.65), and depend on 𝑅𝑠 𝑓𝑓→1, Rgap, 𝑝𝑓𝑓→1 and 𝑧𝑓𝑓→1. From section 2.1.5.2.1, it is known 57 that 𝑝𝑓𝑓→0 ≈ 𝑝𝑓𝑓→1, that is an approximation which is assumed to be able to estimate the other parameters when 𝑓𝑓 → 1. 1 1 0 1 2 ff ct ff ff dl RpC  → → → = (2.64) 1 1 1 1 2 ( ) ff dl ff ff s gap Cz R R  → → → =+ (2.65) Now, the estimation can be performed. The problem is to obtain the value of three parameters from two equations (2.61). As there are more parameters than equations, the choice is to minimize a cost function (CF) to obtain the value of the parameters. The CF must fullfil the condition in such way that real and imaginary parts of (2.60) must be zero. The best way found to meet the condition is to use the module of the complex number ℎ(𝑎𝑜𝑠𝑐, 𝑓 𝑜𝑠𝑐). Due to ℎ(𝑎𝑜𝑠𝑐, 𝑓 𝑜𝑠𝑐)= ℎ1(𝑎𝑜𝑠𝑐, 𝑓 𝑜𝑠𝑐)+ 𝑗 · ℎ2(𝑎𝑜𝑠𝑐, 𝑓 𝑜𝑠𝑐), CF can be defined as: 22 12 ( , ) ( , ) ( , ) osc osc osc osc osc osc h a f h a f h a f+ (2.66) The fact of ℎ1(𝑎𝑜𝑠𝑐, 𝑓 𝑜𝑠𝑐) and ℎ2(𝑎𝑜𝑠𝑐, 𝑓 𝑜𝑠𝑐) are squared assures that terms ℎ1(𝑎𝑜𝑠𝑐, 𝑓 𝑜𝑠𝑐) and ℎ2(𝑎𝑜𝑠𝑐, 𝑓 𝑜𝑠𝑐) cannot compensate each other. Data to fit the model are obtained from a one-week experiment of AA8 cell line with initial number of cells of 2500 cells for wells 1 and 3, 5000 cells for wells 4 and 5, and 10000 cells for wells 7 and 8. Oscillation parameters (𝑓𝐶𝐸 and 𝑎𝐶𝐸) when 𝑓𝑓 → 0 and when 𝑓𝑓 → 1 are summarized in Table 2.1. Notice that 𝑓 𝐶𝐸 𝑓𝑓→0 is the frequency of the signal VCE when 𝑓𝑓 → 0, 𝑓 𝐶𝐸 𝑓𝑓→1 is the frequency of the signal VCE when 𝑓𝑓 → 1, 𝑎𝐶𝐸 𝑓𝑓→0 is the amplitude of the signal VCE when 𝑓𝑓 → 0 and 𝑎𝐶𝐸 𝑓𝑓→1 is the amplitude of the signal VCE when 𝑓𝑓 → 1. Amplitude 𝑎𝐶𝐸 is different from 𝑎𝑜𝑠𝑐 (amplitude of the comparator input voltage) used in DF. The relationship between 𝑎𝑜𝑠𝑐 and 𝑎𝐶𝐸 is 𝑎𝑜𝑠𝑐 = 𝑎𝐶𝐸𝐻𝐾𝐾𝐻𝑃𝐹𝐾𝐿𝑃𝐹 defined in section 2.1.4.2.1. As was explained, there are two steps per well to characterize it. In the first step, using the oscillation parameters for 𝑓𝑓 → 0 (𝑓 𝑜𝑠𝑐 𝑓𝑓→0 and 𝑎𝑜𝑠𝑐 𝑓𝑓→0) and replacing them on the Table 2.1 Oscillation frequency (𝑓𝐶𝐸) and amplitude (𝑎𝐶𝐸) at the beginning of the experiment 𝑓𝑓 → 0 and at the end 𝑓𝑓 → 1 for AA8 cell line. Well 0 [ ] ff CE f Hz → 1 [ ] ff CE f Hz → 0 [ ] ff CE a mV → 1 [ ] ff CE a mV → 1 783.43 867 9.7 13.94 3 802.14 925.43 7.46 15.19 4 821.73 932.12 6.19 13.39 5 780.75 872.2 10.35 14.97 7 822.15 938.3 7.02 16.19 8 843.86 934.37 6.57 14.19 64 ' ' [ ] dlr dl er re dlr dl dl dl e e C C A A C C F C C A A  = → =  (2.72) 4 [ ] 4 e sr e sr s s re re A RA RR RA A     = → =   (2.73) where Cdlr is the interfacial capacitance of the reference electrode model and Rsr is the spreading resistance of the reference electrode model. Note that V0 (equation (2.8)) is different from Cdl (near to 10 mV) and Cdlr (0 V), consequently relationship between Cdl and Cdlr depends on electrode area and the gap resistor (Fig. 2.4), Rgap, is considered equal for microelectrode model and reference electrode model. Since the 𝐴𝑒𝐴𝑟𝑒 ⁄ factor is repeated, the relationship between Ae and Are is defined as: e e re A kA = (2.74) 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 =  = =  (2.75) The reference electrode model is defined as 𝑍𝐶𝐸 𝑟𝑒 (𝑠)= 𝐹(𝑅𝑐𝑡𝑟, 𝐶𝑑𝑙𝑟, 𝑅𝑠𝑟, 𝑅𝑔𝑎𝑝,𝑓𝑓, 𝑠), then TF of the REW model with ten microelectrodes and one reference electrode is: () ( ) ( ) 10 re CE well CE Zs Z s Z s=+ (2.76) where 𝑍𝑤𝑒𝑙𝑙(𝑠) is the impedance of the CE system considering the real distribution of the electrodes. Using the relationships between parameters of the microelectrodes model and reference electrodes model described in (2.75), 𝑍𝑤𝑒𝑙𝑙(𝑠) is a Laplace function that depends on: ( ) ( , , , , , , ) well ct dl s gap e Z s F R C R R ff k s= (2.77) 𝑍𝑤𝑒𝑙𝑙(𝑠) have three poles and three zeros. The expressions of the zeros are too large to present here. In the same way which is done in section 2.1.5.2 poles, zeros and gain for 𝑓𝑓 → 0 and 𝑓𝑓 → 1 must be estimated to fit the model. 00 (1 10 ) (1 10 ) 1 (1 10 ) ct e s e ff ff dl ct dl ct s e R k R k pz CR C R R k →→  +  +  +  − = = −  +  (2.78) 11 (1 10 ) (1 10 ) 11 1 ( (1 10 ) 11 ) ct e s e gap ff ff dl ct dl ct s e gap R k R k R pz CR C R R k R →→  +  +  +  +  − = = −   +  +  (2.79) 0 0 11 10 10 11 11 10 10 10 ff well ct e s e gap ff well ct e s e k R k R k R k R k R k → →     =  + +  +             =  + +  + +          (2.80) 65 In (2.78) are presented the pole and the zero when 𝑓𝑓 → 0, and in (2.79) are presented the pole and the zero when 𝑓𝑓 → 1. These expressions cannot be used to estimate the 𝑍𝑤𝑒𝑙𝑙(𝑠) model parameters, since the approximation done in (2.59) is less accurate for zeros in (2.78) and (2.78). The best choice is to get the 𝑅𝑐𝑡 𝑓𝑓→0 and 𝐶𝑑𝑙 𝑓𝑓→0 expressions from (2.78) and the 𝑅𝑐𝑡 𝑓𝑓→1 and 𝐶𝑑𝑙 𝑓𝑓→1 expressions from (2.79). Working with the equations of poles and zeros, the following expressions are obtained: 00 0 0 (1 10 ) ( ) (1 10 ) s e ff ff ff ct ff e R k p z Rpk →→ → →  +   − =−  +  (2.81) 0 00 1 10 (1 10 ) ( ) ff e dl s e ff ff k CR k p z → →→ + =−  +   − (2.82) 00 1 0 ( (1 10 ) 11 ) ( ) (1 10 ) s e gap ff ff ff ct ff e R k R p z Rpk →→ → →  +  +   − =−  +  (2.83) 1 00 1 10 ( (1 10 ) 11 ) ( ) ff e dl s e gap ff ff k CR k R p z → →→ + =−  +  +   − (2.84) Then, the values of Rct, Cdl, poles and zeros when 𝑓𝑓 → 0 and 𝑓𝑓 → 1 can be obtained. The values of 𝑓 𝐶𝐸 𝑓𝑓→0, 𝑓 𝐶𝐸 𝑓𝑓→1, 𝑎𝐶𝐸 𝑓𝑓→0 and 𝑎𝐶𝐸 𝑓𝑓→1 used in the step one and the step two of the fitting process are the same of the Table 2.1. In Table 2.7 and Table 2.8 the values of poles, zeros and model parameters are summarized. Note that, as with the SEW model, 𝑝𝑓𝑓→0 ≈ 𝑝𝑓𝑓→1 have been considered. 𝑅𝑐𝑡 𝑓𝑓→0 and 𝑅𝑐𝑡 𝑓𝑓→1 values are still an order of magnitude below 13.67 MΩ expected from (2.6). Values of 𝐶𝑑𝑙 𝑓𝑓→0 and 𝐶𝑑𝑙 𝑓𝑓→1 are very close to the expected theoretical value 3.29 nF in (2.11). Besides that, Rs values are different from the expected value of 1 kHz in (2.13). The value of Rgap is studied below. fval of both fitting steps for each well is shown in Table 2.9. Notice that, compared with fval for SEW model (Table 2.2), 2 less orders of magnitude have been achieved. As for SEW model, cell concentration is estimated using ff for each moment and well. ff estimation is perform in the same way of SEW model, but in that case, the REW model Table 2.7 Model parameters at the beginning of the experiment 𝑓𝑓 → 0 for AA8 cell line (REW model with 10 microelectrodes and one reference electrode). Well 0 [ ] ff ct RM → 0 [ ] ff dl C nF → 0 [ ] ff s Rk → 0 [ ] ff p Hz → 0 [ ] ff z kHz → 1 1.17 2.74 4.17 49.67 9.59 3 0.71 3.75 3.23 60.17 9.04 4 0.52 4.8 2.95 64.24 7.75 5 1.41 2.54 4.71 44.2 9.12 7 0.63 4.24 3.44 59.78 7.52 8 0.36 4.85 2.58 91.65 8.79 66 of the wells is used. Limits for CF minimization are the same of Table 2.5, but the values of 𝑅𝑐𝑡 𝑓𝑓→0, 𝑅𝑐𝑡 𝑓𝑓→1, 𝐶𝑑𝑙 𝑓𝑓→0 and 𝐶𝑑𝑙 𝑓𝑓→1 used for each well are the values of Table 2.7 and Table 2.8. Table 2.10 present the values for Rsi and ΔRs for the REW model. In Fig. 2.25 shows the fill-factor evolution for each well using the parameters obtained in fitting process. As expected, wells 1 and 3 are the last ones to reach the maximum fillfactor, wells 4 and 5 are the second ones and the first to reach the maximum are wells 7 and 8. This is due to the initial number of cells on each well. The more the initial number of cells, the lower the time to reach the maximum fill-factor. In the same way for SEW model, results of ff can be used to calculate cell concentration using (2.70). It is interesting to compare SEW model result with REW model results and with cell concentration of traditional counting method. Fig. 2.26 shows Table 2.8 Model parameters at the end of the experiment 𝑓𝑓 → 1 for AA8 cell line (REW model with 10 microelectrodes and one reference electrode). Well 1 [ ] ff ct RM → 1 [ ] ff dl C nF → 1 [ ] ff s Rk → [ ] gap R 1 [ ] ff z kHz → 1 1.28 2.5 6.08 346.6 5.12 3 0.88 3.02 7.41 383.5 3.63 4 0.69 3.61 1.6 1021.1 3.36 5 1.51 2.38 7.79 258.5 4.8 7 0.87 3.07 5.01 841.6 3.24 8 0.51 3.43 1.2 1017.4 3.73 Table 2.9 Function value (𝑓𝑣𝑎𝑙) for each well in the two steps of model parameters estimation, for 𝑓𝑓 → 0 and 𝑓𝑓 → 1. Well 1 3 4 5 7 8 0ff val f→ [x 10-5] 0.089 0.38 0.99 0.161 0.16 0.69 1ff val f→ [x 10-5] 0.26 0.068 1.53 0.18 1.61 0.38 Table 2.10 Parameters to define the evolution of Rs. Well [ ] si Rk [ ] s Rk 1 4.17 1.91 3 3.23 4.18 4 2.95 -1.35 5 4.71 3.09 7 3.44 1.57 8 2.58 -1.38 67 this comparison. Cell concentration evolution of REW model seems to be similar to cell concentration evolution of traditional counting method, except for the wells with an initial number of cells of 10000, where it differs slightly. In general, REW model seems to show better the behavior (or cell concentration evolution) of cell concentration (or fillfactor) for low values of ff. Like the SEW model, the REW model meets the error margins obtained from traditional optical counting. In this subsection, it has been demonstrated that linear models can be used to know the behavior of a CC assay even if it has a clearly nonlinear behavior. Therefore, it can be concluded that it is possible to obtain the cell concentration of a CC assay using an oscillator circuit (OBT) as a measurement circuit. The errors found are relatively small, and the trend of the traditional curves and those using the OBT circuit are similar. a) b) c) Fig. 2.25 Fill-factor evolution for each well of the experiment using the REW model. a) Initial number of cells: 2500 cells (well 1 on the left and well 3 on the right). b) Initial number of cells: 5000 cells (well 4 on the left and well 5 on the right). c) Initial number of cells: 10000 cells (well 7 on the left and well 8 on the right). 68 050 100 150 Time [h] 0 2 4 6 8 10 12 14 16 18 Cell concentration [cells/cm2] 104 Nini = 2500 cells CF minimization Real-electrode well CF minimization Single-electrode well Traditional method 050 100 150 Time [h] 0 0.5 1 1.5 2 2.5 Cell concentration [cells/cm2] 105 Nini = 5000 cells CF minimization Real-electrode well CF minimization Single-electrode well Traditional method 050 100 150 Time [h] 0 0.5 1 1.5 2 2.5 3 Cell concentration [cells/cm2] 105 Nini = 10000 cells CF minimization Real-electrode well CF minimization Single-electrode well Traditional method Fig. 2.26 Cell concentration comparison between CF minimizing method using SEW model (blue) and REW model (red), and traditional counting method (yellow). 69 2.1.5.4 Fractional order model Integer-Order (IO) models, such as those used in the previous subsections, do not completely fit the behavior of a real system. In the case of the electrical behavior of a BS, their BDs confirm that a (relatively simple) IO model cannot capture all the dynamics of the response to an electrical stimulus. Although in the previous subsections the evolution of the growth of a CC assay has been successfully obtained using IO models (with a certain degree of complexity), it is interesting to apply Fractional Order (FO) models in the process of extracting ff from osc f and osc a . FO models can achieve a better fit of the CE system parameters, because they are non-linear and are closer to the real behavior. For example, in [35] a two electrode model, which contains a FO Constant Phase element (CPE) is used to is used to perform the characterization of electrodes for bioimpedance measurement of animal tissue. Authors consider that the electrodes have resistive and capacitive properties of FO (the CPE element), and that the tissue is of the resistive type. FO elements can also be used to characterize vegetal tissues, as is done in [34]. FO models are based on the premise that the order of a differential operator can be non-integer. To start learning about FO elements and their frequency behavior, the manual [64] is used as a reference document. The authors define the FO Laplace operator as follows, 0< <2 n nn s   = (2.85) where αn is the order of the differential Laplace operator s and n is an integer such that n = {1, 2, 3…}. Capacitors of the CE model (Fig. 2.5) are replaced by CPE described with FO operators, 12 12 12 11 dl dl CC dl dl XX s C s C  == (2.86) where α1 is the order of s in reactance 𝑋𝐶𝑑𝑙1 and α2 is the order of s in reactance 𝑋𝐶𝑑𝑙1. The transfer function of the 10 electrodes FO CE model is very complex, so it is expressed as: 12 ( ) ( , , , , , , , , ) well ct dl s gap e Z s F R C R R ff k s  = (2.87) From first tests, where a CF of the FO CE model was minimized to know the possible values of α1 and α2, a conclusion is reached: in this FO CE model, Rs does not change with ff (i.e., ΔRs is removed from the model). The first tests have reported values of ΔRs in the order of 10-3, which does not affect the model because Rs has an order of 103. The explanation begins by derived the FO models of the 10 electrodes system when 𝑓𝑓 → 0 and 𝑓𝑓 → 1, in order to estimate the parameters in a simple way. The transfer functions obtained are as follows, 70 0 1 0 00 1 0 1 0 () ff ff ff ff well well ff ff pz Zk zp   →→ →→ →→ + =+ (2.88) 1 2 1 11 2 1 2 1 () ff ff ff ff well well ff ff pz Zk zp   →→ →→ →→ + =+ (2.89) where poles, zeros and gains have been previously defined in equations (2.78) to (2.80). Note that, when 𝑓𝑓 → 0, α2 has no influence in the system behavior, and, when 𝑓𝑓 → 1, α1 has any influence in the system behavior. In an FO model, α changes the slope of poles and zeros, such that its magnitude curve has a slope of ± α·20 dB/dec. Furthermore, the change in phase, for single poles and zeros, is ± α·π/2 º. The experimental BDs (Fig. 2.19) of the system indicate that the system has a non-linear behavior, because the slope of the magnitude curve and the changes in the phase curve do not match the general characteristics of an IO model. Considering that the fractional exponent α can modify these indicators outside the standard of an IO model, and knowing that Rs defines zeros and gains in the IO CE model, the conclusion is that small changes between α1 and α2 (transition from 𝑓𝑓 → 0 to 𝑓𝑓 → 1) produce big changes in the FO model of the system, and, therefore, ΔRs becomes redundant and can be removed from the FO model. The parameter estimation process is similar to the performed for the 1 electrode and 10 electrodes IO models. The BSC is used to ensure the oscillation conditions, obtaining a complex number (2.60), which must be equal to 0. Since s = j·ω, in FO models complex number j is replaced by jα. Then, the FO complex number must be translated to complex number by using Moivre's formula: cos sin 22 jj       =  +           (2.90) Due to the high complexity of the model (cross products of α1 and α2), the model is implemented in a different way than the previous sections. IO models have been implemented (sections 2.1.5.2 and 2.1.5.3) completely to obtain its parameters, while FO models are, in this case, implemented in a transitional mode. This means that, considering that the pole is constant for any value of ff, a transition from 𝑧𝑓𝑓→0 to 𝑧𝑓𝑓→1 and from 𝑘𝑤𝑒𝑙𝑙 𝑓𝑓→0 to 𝑘𝑤𝑒𝑙𝑙 𝑓𝑓→1 is implemented using ff. In other words, ff is used to change from 𝑍𝑤𝑒𝑙𝑙 𝑓𝑓→0(𝜆1) to 𝑍𝑤𝑒𝑙𝑙 𝑓𝑓→1(𝜆2). Then, the implemented model is: 1, 2 () well well pz Zk zp    + =+ (2.91) where 𝑝 = 𝑝𝑓𝑓→0 = 𝑝𝑓𝑓→1 and, 01 (1 ) ff ff z ff z ff z →→ = −  +  (2.92) 01 (1 ) ff ff well well well k ff k ff k →→ = −  +  (2.93) 12 (1 )ff ff    = −  +  (2.94) 71 When 𝑓𝑓 → 0, by replacing fosc by 𝑓 𝐶𝐸 𝑓𝑓→0 and aosc by 𝑎𝐶𝐸 𝑓𝑓→0 (Table 2.1) in (2.91), 𝑅𝑐𝑡, 𝐶𝑑𝑙, 𝑅𝑠, 𝑝, 𝑧𝑓𝑓→0 and α1 can be obtained by minimizing the CF. The CF for FO model is the same of (2.66), except that during its calculation the CE model used is the FO real electrode-well (REW) model (2.91). In the same way, when 𝑓𝑓 → 1, by replacing fosc by 𝑓 𝐶𝐸 𝑓𝑓→1 and aosc by 𝑎𝐶𝐸 𝑓𝑓→1 (Table 2.1) in (2.91), 𝑅𝑔𝑎𝑝, 𝑧𝑓𝑓→1 and α2 can be obtained by minimizing the CF (2.66). As in the previous sections, data obtained from the AA8 cell line is used. In this case, to minimize the CF, the Matlab® function fminsearch [65] is used. It performs an unbounded search of the values that make the CF minimum, starting from the initial point chosen by the user. The value of the function, fval, is again used as a measurement of the minimization quality. Table 2.11 shows the fval values obtained for each well and ff limit. Notice that fminsearch finds parameter values for which fval is very small. The results of minimizing the CF of (2.66) are shown in next tables. Table 2.12 displays the CE model parameters when 𝑓𝑓 → 0 and Table 2.13 presents the CE model Table 2.11 Function value (𝑓𝑣𝑎𝑙) for each well in the two steps of FO model parameters estimation, for 𝑓𝑓 → 0 and 𝑓𝑓 → 1. Well 1 3 4 5 7 8 0ff val f→ [x 10-15] 0.4 0.18 0.72 0.04 0.38 0.26 1ff val f→ 0 0 0 0 0 0 Table 2.12 Model parameters at the beginning of the experiment 𝑓𝑓 → 0 for AA8 cell line (FO REW model with 10 microelectrodes and one reference electrode). Well [ ] ct RM [ ] dl C nF [ ] s Rk [ ]p Hz 0 [ ] ff z kHz → 1  1 0.43 1.85 0.97 200 61.1 1.05 3 0.29 2.72 0.29 199.1 137.65 1.04 4 0.47 4.1 2.39 83.05 11.2 1.02 5 0.69 2.42 0.56 94.7 80.63 1 7 0.24 3.03 0.36 222.98 98.89 1.04 8 0.84 7.2 1.21 26.24 12.52 0.95 Table 2.13 Model parameters at the end of the experiment 𝑓𝑓 → 1 for AA8 cell line (FO REW model with 10 microelectrodes and one reference electrode). Well [ ] gap R 1 [ ] ff z kHz → 2  1 501.13 12.89 1.03 3 699.81 7.58 1 4 520.64 4.37 0.98 5 686.71 8.12 0.99 7 355.63 12.32 0.99 8 412.02 3.6 0.91 72 parameters when 𝑓𝑓 → 1. The transition from α1 to α2 models the effect of changes in the zero and gain. Thus, Rct, Cdl and Rs do not change during the CC assay growth, they are obtained at 𝑓𝑓 → 0 and are constant. Fig. 2.27 shows ff evolution for each well using the parameters obtained in fitting process. Results are similar to the ones in subsections 2.1.5.2 and 2.1.5.3. In this case, it seems that ff tends to higher values in the first days of simulation. Even so, the trend and the time at which ff reaches value 1 is similar to the previous subsections. As in previous sections, ff is used to calculate cell concentration using (2.70). Cell concentration results of FO REW model are once again compared with SEW model results and REW model results and with cell concentration of traditional counting method. Fig. 2.28 shows this comparison. Cell concentration evolution of FO REW model seems to be similar to cell concentration evolution of the previous models. The results prove that FO models can be used to track the growth of CC assays, but with a sufficiently high initial ff or cell concentration. For Nini values of 2500 and 5000 cells, it has significant errors in the first hours and/or days of the experiment. All three models (SEW, REW and FO REW) predict correctly the growth trend of the CC assay. However, the FO REW model, for low values of Nini, fits worse than the IO models, which may be due to the high sensitivity of α1 and α2. When a FO is minimized, several local minima can be found, making it difficult for the algorithm to know which is the correct value. On the other hand, FO model fits better the BD of the real system, a) b) c) Fig. 2.27 Fill-factor evolution for each well of the experiment using the FO REW model. a) Initial number of cells: 2500 cells (well 1 on the left and well 3 on the right). b) Initial number of cells: 5000 cells (well 4 on the left and well 5 on the right). c) Initial number of cells: 10000 cells (well 7 on the left and well 8 on the right). 73 050 100 150 Time [h] 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Cell concentration [cells/cm2] 105 Nini = 2500 cells CF minimization Real-electrode well (Fractional Order) CF minimization Real-electrode well CF minimization Single-electrode well Traditional method 050 100 150 Time [h] 0 0.5 1 1.5 2 2.5 Cell concentration [cells/cm2] 105 Nini = 5000 cells CF minimization Real-electrode well (Fractional Order) CF minimization Real-electrode well CF minimization Single-electrode well Traditional method 050 100 150 Time [h] 0 0.5 1 1.5 2 2.5 3 Cell concentration [cells/cm2] 105 Nini = 10000 cells CF minimization Real-electrode well (Fractional Order) CF minimization Real-electrode well CF minimization Single-electrode well Traditional method Fig. 2.28 Cell concentration comparison between CF minimizing method using FO REW model (blue), REW model (red) and SEW model (yellow), and traditional counting method (purple). 80 In Table 2.21 to Table 2.24 are shown the results of fitting the N2aAPP cell line with the REW model. The results are similar to those obtained for N2a, with a few more dispersion in some parameters. Table 2.21 Model parameters at the beginning of the experiment 𝑓𝑓 → 0 for N2aAPP cell line (REW model with 10 microelectrodes and one reference electrode). Well 0 [ ] ff ct RM → 0 [ ] ff dl C nF → 0 [ ] ff s Rk → 0 [ ] ff p Hz → 0 [ ] ff z kHz → 2 0.61 3.20 2.38 81.24 14.35 3 0.43 3.07 1.48 119.55 24.09 4 0.35 3.99 1.67 115.30 16.42 6 0.47 3.65 1.58 93.19 18.86 7 0.43 3.75 1.62 97.81 17.94 8 0.32 4.33 1.91 114.08 13.26 Table 2.22 Model parameters at the end of the experiment 𝑓𝑓 → 1 for N2aAPP cell line (REW model with 10 microelectrodes and one reference electrode). Well 1 [ ] ff ct RM → 1 [ ] ff dl C nF → 1 [ ] ff s Rk → [ ] gap R 1 [ ] ff z kHz → 2 0.66 2.96 1.75 417.65 7.67 3 0.48 2.78 0.49 628.83 7.76 4 0.39 3.57 2.23 384.32 6.17 6 0.39 4.39 4.72 69.89 4.91 7 0.47 3.45 2.66 162.74 8.35 8 0.29 4.71 3.14 415.53 3.87 Table 2.23 Parameters to define the evolution of Rs. Well [ ] si Rk [ ] s Rk 1 2.38 -0.63 3 1.48 -0.99 4 1.67 0.56 5 1.58 3.14 7 1.62 1.04 8 1.91 1.22 Table 2.24 Function value (𝑓𝑣𝑎𝑙) for each well in the two steps of model parameters estimation of N2a, for 𝑓𝑓 → 0 and 𝑓𝑓 → 1. Well 2 3 4 6 7 8 0ff val f→ [x 10-15] 0 0.11 0.56 0 0 0 1ff val f→ [x 10-5] 0.35 0.31 1.44 0.15 0.28 0.51 81 In the same way, the results of the FO REW model fit for the N2aAPP cell line are similar to those obtained for the N2a line, and are shown in Table 2.25 to Table 2.27. All the data will be discussed in the following section in which an explanation of the dispersions will be proposed by relating them to various parameters of cell morphology and the oscillation waveforms (𝑓𝐶𝐸 and 𝑎𝐶𝐸). Fig. 2.31 and Fig. 2.32 shows the cell concentration comparison of N2aAPP cell line. Cell concentration obtained by using parameters values of IO and FO models, is compared with the cell concentration obtained from the traditional optical cell count. As for the N2a line, the cell growth time of the traditional methods is longer than that obtained from the prototype data (Fig. 2.31). Therefore, in Fig. 2.32, the time of cell concentration obtained by traditional optical counting is advanced in order to compare the evolution of the concentrations. For low concentrations, the model also does not Table 2.25 Model parameters at the beginning of the experiment 𝑓𝑓 → 0 for N2aAPP cell line (FO REW model with 10 microelectrodes and one reference electrode). Well 0 [ ] ff ct RM → 0 [ ] ff dl C nF → [ ] s Rk 0 [ ] ff p Hz → 0 [ ] ff z kHz → 1  2 0.90 6.21 1.36 28.49 12.90 0.95 3 1.44 6.65 1.57 16.63 10.43 0.94 4 1.13 6.46 1.24 21.80 13.55 0.93 6 0.37 4.97 1.81 87.05 12.15 0.99 7 0.36 5.19 1.55 85.83 13.59 0.99 8 0.26 5.41 1.22 111.76 16.60 0.99 Table 2.26 Model parameters at the end of the experiment 𝑓𝑓 → 1 for N2a cell line (FO REW model with 10 microelectrodes and one reference electrode). Well 1 [ ] ff ct RM → 1 [ ] ff dl C nF → [ ] gap Rk 1 [ ] ff z kHz → 2  2 3.75 1.49 1.46 6.07 1.13 3 2.14 4.47 0.65 3.84 0.99 4 3.59 2.03 1.47 4.45 1.08 6 0.63 2.92 0.92 4.43 1.06 7 0.84 2.22 1.11 5.14 1.09 8 0.54 2.64 0.92 5.29 1.09 Table 2.27 Function value (𝑓𝑣𝑎𝑙) for each well in the two steps of FO model parameters estimation, for 𝑓𝑓 → 0 and 𝑓𝑓 → 1. Well 1 3 4 5 7 8 0ff val f→ [x 10-15] 0.34 0.32 0.11 0.23 0.23 0.22 1ff val f→ [x 10-5] 1.72 0.35 1.94 0.42 0.44 1.55 82 properly reflect the real behavior of the CC assay. The trend of the curves is similar for medium and high concentration values, but like the N2a line, the fit is not sufficiently accurate, as it exceeds the error margins at many points. The data from the OBT oscillation connected to the CC assays of the N2a and N2aAPP cell lines have provided evidence that the method of model fitting and subsequent cell concentration estimation works. Notice that, both cell lines present a worse match than the AA8 line, and therefore this point should be analyzed. 83 050 100 150 Time [h] 0 1 2 3 4 5 6 7 8 9 Cell concentration [cells/cm2] 105 Nini = 2500 cells CF minimization Real-electrode well (Fractional Order) CF minimization Real-electrode well Traditional method 050 100 150 Time [h] 0 1 2 3 4 5 6 7 8 9 Cell concentration [cells/cm2] 105 Nini = 5000 cells CF minimization Real-electrode well (Fractional Order) CF minimization Real-electrode well Traditional method 050 100 150 Time [h] 0 1 2 3 4 5 6 7 8 9 10 Cell concentration [cells/cm2] 105 Nini = 10000 cells CF minimization Real-electrode well (Fractional Order) CF minimization Real-electrode well Traditional method Fig. 2.31 N2aAPP cell line cell concentration comparison between CF minimizing method using FO REW model (blue) and REW model (red), and traditional counting method (yellow ). 84 010 20 30 40 50 60 70 80 90 100 Time [h] 0 1 2 3 4 5 6 7 8 Cell concentration [cells/cm2] 105 Nini = 2500 cells CF minimization Real-electrode well (Fractional Order) CF minimization Real-electrode well Traditional method 010 20 30 40 50 60 70 80 90 100 Time [h] 0 1 2 3 4 5 6 7 8 9 Cell concentration [cells/cm2] 105 Nini = 5000 cells CF minimization Real-electrode well (Fractional Order) CF minimization Real-electrode well Traditional method 010 20 30 40 50 60 70 80 90 100 Time [h] 0 1 2 3 4 5 6 7 8 9 10 Cell concentration [cells/cm2] 105 Nini = 10000 cells CF minimization Real-electrode well (Fractional Order) CF minimization Real-electrode well Traditional method Fig. 2.32 N2aAPP cell line cell concentration comparison between CF minimizing method using FO REW model (blue) and REW model (red), and traditional counting method (yellow). The time vector has been advanced by 48 hours to match the lines. 85 2.1.5.5.2 Results comparison between cellular lines The previous sections show that the results obtained, using cell concentration as a metric, are better for the AA8 cell line than for N2a and N2aAPP. In this section, an explanation for this will be discussed, as it is interesting to know the reason for future experiments. The first factor to consider is the frequency and amplitude of oscillation. The metrics to be used are 𝑓𝐶𝐸 and 𝑎𝐶𝐸. Concretely, it is observed that the difference between 𝑓𝐶𝐸 and 𝑎𝐶𝐸 between 𝑓𝑓 → 0 and 𝑓𝑓 → 1 is greater in the AA8 line than in the other cell lines fitted. If 𝑓 𝐶𝐸 𝑓𝑓→0 and 𝑎𝐶𝐸 𝑓𝑓→0are defined as the frequency and amplitude when 𝑓𝑓 → 0, and 𝑓 𝐶𝐸 𝑓𝑓→1 and 𝑎𝐶𝐸 𝑓𝑓→1 are defined as the frequency and amplitude when 𝑓𝑓 → 1, the variations of 𝑓𝐶𝐸 and 𝑎𝐶𝐸 can be defined as 10 10 ff ff CE CE CE ff ff CE CE CE f f f a a a →→ →→  = −  = − (2.95) where 𝛥𝑓𝐶𝐸 and 𝛥𝑎𝐶𝐸 are the difference in frequency and amplitude, respectively, between the time when 𝑓𝑓 → 1 and 𝑓𝑓 → 0. Table 2.28 and Table 2.29 shows the values of 𝛥𝑓𝐶𝐸 and 𝛥𝑎𝐶𝐸, respectively, of the three cell lines. The tables reveal a greater difference in amplitude variation between the AA8 cell line and the other two than in frequency. To find a relationship with the BI, the magnitude and phase variation of the BI from 𝑓𝑓 → 0 to 𝑓𝑓 → 1 is calculated. As for frequency and amplitude, the magnitude and phase variation can be defined as 10 10 ff ff CE CE CE ff ff CE CE CE Mag Mag Mag Ph Ph Ph →→ →→  = −  = − (2.96) Table 2.28 𝛥𝑓𝐶𝐸 values of AA8, N2a and N2aAPP cell lines [Hz] Nini 2500 5000 10000 AA8 83.57 90.51 91.45 110.39 116.15 123.28 N2aAPP 52.15 70.26 53.37 104.92 60.95 111.64 N2a 56.15 69.99 61.22 75.14 67.22 103.40 Table 2.29 𝛥𝑎𝐶𝐸 values of AA8, N2a and N2aAPP cell lines [mV] Nini 2500 5000 10000 AA8 4.23 4.61 7.20 7.63 7.73 9.17 N2aAPP 2.04 3.65 2.27 3.66 2.41 3.80 N2a 0.96 1.92 1.10 1.98 1.27 2.97 86 where 𝛥𝑀𝑎𝑔𝐶𝐸 and 𝛥𝑃ℎ𝐶𝐸 are the difference in magnitude and phase, respectively, between the time when 𝑓𝑓 → 0 and 𝑓𝑓 → 1. Table 2.30 and Table 2.31 show the values of 𝛥𝑀𝑎𝑔𝐶𝐸 and 𝛥𝑃ℎ𝐶𝐸 sorted by cellular line and Nini. 𝛥𝑀𝑎𝑔𝐶𝐸 has different values depending on the cell line, whereas 𝛥𝑃ℎ𝐶𝐸 takes close values for different cell lines. Due to the highest differences in the oscillation parameters between cell lines, 𝛥𝑀𝑎𝑔𝐶𝐸 and 𝛥𝑃ℎ𝐶𝐸 are presented versus to the 𝛥𝑎𝐶𝐸 for each cell line (note that the data are sorted in ascending order of 𝛥𝑎𝐶𝐸). In this way, Fig. 2.33 shows these relationships, where data of Table 2.30 and Table 2.31 is displayed versus data of Table 2.29. In the figures there are three clearly defined zones, one for each cell line. This clearly suggests that the cells (size, shape, composition, etc.) influence the accuracy of the fitting method. A useful parameter, known from experiments by traditional methods (optical counting), is the mean area of the CC assay cells. The area of the cells affects 2 main parameters: Rgap and Rs. The Table 2.32 shows the mean area of each cell line and the ranges of Rgap and Rs obtained for each cell line for the REW model. In [55] the author of the present thesis shown how both parameters influence 𝑎𝐶𝐸 and 𝑓𝐶𝐸 when 𝑓𝑓 → 0 to 𝑓𝑓 → 1. On the other hand, looking at Fig. 2.4, it is concluded that Rgap and Rs (Rgap in a stronger degree) depend on the size and morphology of the cells. So, these facts indicate that the cell area influences 𝛥𝑎𝐶𝐸 and 𝛥𝑓𝐶𝐸, which results in a loss of accuracy in the model fitting method when the cell area is low. Table 2.30 𝛥𝑀𝑎𝑔𝐶𝐸 values of AA8, N2a and N2aAPP cell lines [V/V] Nini 2500 5000 10000 AA8 30.16 654.77 993.03 -209.67 1288.03 1381.04 N2aAPP 166.03 -1374.60 329.02 180.28 271.68 -621.72 N2a -577.87 -372.62 -650.96 -384.35 -643.13 509.08 Table 2.31 𝛥𝑃ℎ𝐶𝐸 values of AA8, N2a and N2aAPP cell lines [°] Nini 2500 5000 10000 AA8 4.57 8.72 8.17 5.10 9.06 7.22 N2aAPP 2.63 7.42 3.72 2.87 4.68 9.04 N2a 3.26 7.22 3.59 5.20 4.58 4.83 Table 2.32 Values of Acell, and ranges of Rgap and Rs of the cell lines under study for REW model. Cell line Acell [μm2] Rgap [Ω] 0 [ ] ff s Rk → AA8 552.5 258.5 – 1021.1 2.58 – 4.71 N2aAPP 118 69.9 – 628.8 1.48 – 2.38 N2a 184 43.1 – 512.7 1.39 – 1.72 87 Table 2.33 shows the same comparison of Table 2.32, but for FO REW model. As can be seen in the cell concentration figures (Fig. 2.28 for AA8, Fig. 2.30 for N2a and Fig. 2.32 for N2aAPP), FO model fits better to the results obtained by traditional counting than the IO model. It can be concluded that, for this particular case, the FO model captures better the growth of a CC assay, since the model parameters are closer to the real ones because 0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01 aCE [mV] -1500 -1000 -500 0 500 1000 1500 Relationship between aCE and MagCE AA8 N2a N2aAPP MagCE [V/V] 0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01 aCE [mV] 2 3 4 5 6 7 8 9 10 Relationship between aCE and PhCE AA8 N2a N2aAPP PhCE [°] Fig. 2.33 Relationship between 𝛥𝑎𝐶𝐸 and oscillation parameters. Top figure shows the relationship 𝛥𝑎𝐶𝐸 - 𝛥𝑀𝑎𝑔𝐶𝐸, and bottom figure shows the relationship 𝛥𝑎𝐶𝐸 - 𝛥𝑃ℎ𝐶𝐸. Table 2.33 Values of Acell, and ranges of Rgap and Rs of the cell lines under study for FO REW model Cell line Acell [μm2] Rgap [Ω] [ ] s Rk AA8 552.5 355.6 – 699.8 0.29 – 2.39 N2aAPP 118 493.8 – 2010.6 1.12 – 1.98 N2a 184 653 – 1474.2 1.22 – 1.81 88 they should not compensate for the nonlinearities modeled by the FO Laplace operator exponent (2.85). Comparing the three cell lines, it can be seen that while AA8 cells are easily fitted, the other two lines are not. This is mostly due to the range of the oscillation amplitude. The (2.35) shows the dependence between the oscillation amplitude and the BI of the CE system. The simplest solution to this problem is to tune the oscillator to a different frequency than 1 kHz. The BD of the AA8 cellular line is very useful to know the optimum frequency. Fig. 2.34 shows the BD of well in a AA8 cell line experiment. Magnitude and phase changes with ff, but for the tunned frequency (1 kHz) the variation is much smaller than for higher frequencies (10 kHz). In the absence of experimental BDs of the N2a and N2aAPP lines, the optimal frequency to maximize the range cannot be known, but the BD of Fig. 2.34 suggests that the accuracy of the tuning method improves greatly when the range (which depends on ff) is higher. Therefore, for future experiments the bandpass filter will be tuned to higher frequencies. The second problem faced during the fitting of the N2a and N2aAPP cell lines is the growth rate of the CC assay. While the culture of the AA8 line grows at the same rate as the traditional parallel experiment, the other two lines seem to grow much faster than the traditional experiment performed for the N2a and N2aAPP lines. Taking into account that the conditions of the traditional experiments and using the OBT are the same, that the initial cells have the same origin, and the experiments have been performed in parallel, the only possible explanation is electrostimulation. Since the current through the cells is similar for all cell lines (close to 7 μA), the only difference is the cells themselves. Table 2.32 shows the average cell sizes for each cell line. The hypothesis proposed in this work is that the current increases the growth rate, depending on the size of the cells. Therefore, the AA8 cell line, having a considerably Δff Δff Fig. 2.34 BD of a well in AA8 cell line experiment. The range of the magnitude and phase changes with the oscillation frequency. 89 larger size (on average) than the other two lines, is not affected by the current levels used, while the N2a and N2aAPP lines are accelerated in growth due to the currents used. This is only a hypothesis, which is still under study, whose results will be studied in the future. 2.2 Conclusions This chapter has explained in detail how to model the behavior of the growth of a CC assay using OBT data. The oscillator changes its amplitude and frequency of oscillation as a function of the cell concentration from the cell culture which is connected as a load. Some conclusions and questions arising from the modeling work are detailed below: • It is possible to measure the growth of a CC assay using an oscillator. The BS is connected as a load, which by varying its BI, changes the amplitude and frequency of oscillations. In addition, using the oscillation frequency and amplitude, and the BSC, the electrical model of the BS can be known, while injecting a signal of a given frequency through a CC assay, only its impedance at the frequency of the injected signal is known, so it would be necessary to make a sweep in frequency to obtain the BI model. • The DF method is a very accurate approach to linearize a non-linear element, which allows to use the BSC to obtain the theoretical (mathematically) oscillation frequencies and amplitudes. • A CC assay can be modeled by its BI, which is highly non-linear. Therefore, the theoretical value of a parameter of a CE electrical model should be considered as a nominal value with a certain margin of error. • Experience with the three cell lines used in this work suggests that a tunable BPF is necessary in order to be able to adapt the oscillation frequency to the most suitable zone according to the electrical behavior of the BS. • The SEW model allows easy fitting of a CC assay. On the other hand, it does not work sufficiently robust for low Nini, or in case the oscillation parameters have low ranges (e.g., the N2a and N2aAPP lines). • The experimental BD is very useful to configure the measurement circuit (OBT), since it provides the optimum oscillation frequency for each BS. • Choosing an appropriate CF is important, as the minimization of CF must meet all the objectives required. For that reason, the module of an imaginary number is perfect as a CF for the fitting which use the BSC. 96 34 34 1 2 ( ) ( 1) 10 100 10 10 10 10 0 1 ini ini gap ini ff ini s ini ff j ff j R z R  → =−  =  = = = (3.3) The steps described above are applied for each value of j for each well of each cell line, performing simulations that do not consider the future values of 𝑓𝐶𝐸 and 𝑎𝐶𝐸, using the REW and FO REW models. The following subsection shows the results of this simulation procedure of the real CC assay experiments, which predicts the cell concentration in RT. 3.1.2 Real-time simulation results The RT simulation method is finally designed to be implemented in a prototype to report the ff and cell concentration values after each measurement. As the required sampling time is 1 h, the time to perform all the necessary mathematical operations is not a critical point. This means that the computation time of the initial parameters, and of the ff and other parameters for each measurement, is not a critical problem. Two metrics will be used to determine the level of accuracy of the method. The first is the error in the estimation of ff, and the second is the error in the calculation of the cell concentration. The error in the estimation of the fill factor will be measured with respect to the deviation from its ideal final value. It is considered that for the maximum ff obtained in the simulation, the error is the difference between the maximum value obtained and the maximum value that ff should reach, i.e., a value of 0.99. Thus, it is possible to calculate the real ff curve that should have been obtained for each well using the following equation, ( ) ( ) ( ) 0.99 max i ff sim ff sim ff j k ff j kff = = (3.4) where ffsim is the vector of ff values obtained from RT simulation, ffi is the vector of theoretically real values of ff and kff is the factor applied to ffsim to obtain ffi. The metrics for measuring accuracy are the mean relative error in the fill factor and the mean relative error in cell concentration, both in percent (%).The results for each of the cell lines are presented in the next three subsections. 3.1.2.1 AA8 results This cell line is the one that should obtain better results, considering the Chapter 2, where the best results of modelling CE block were achieved with the AA8 cell line. The 97 following tables show the values of the initial parameters of the simulation in RT, i.e., when the simulator detects that the minimum 𝑓𝐶𝐸 and 𝑎𝐶𝐸 have been reached. Table 3.1 shows the values obtained using the REW model, while Table 3.2 shows the values obtained using the FO REW model. Large differences are observed between some of the parameter values of the two models, for example, between the parameters 𝑅𝑐𝑡 𝑓𝑓→0, 𝑝𝑓𝑓→0 and 𝑧𝑓𝑓→0 of both models. This effect may be due to the fact that by introducing α1, the other parameters do not have to adjust for the nonlinear behavior of the real CE system. If the obtained values for the REW model (Table 3.1) are compared with the values obtained for the same model during the modeling (Table 2.7), it is observed that the new values have less dispersion. The differences during the estimation of the initial parameters are that 𝑓 𝐶𝐸 𝑓𝑓→0 and 𝑎𝐶𝐸 𝑓𝑓→0 are obtained automatically, and that the value of ff parameter is estimated at the same time as the other parameters. The value of ff when 𝑓𝐶𝐸 and 𝑎𝐶𝐸 reach its minimum is named as ffi. The fact that ffi is estimated, gives more freedom to the other parameters to get closer to their correct value. Table 3.1 Model parameters estimated at the beginning of the RT simulation (𝑓𝑓 → 0) for AA8 cell line (REW model with 10 microelectrodes and one reference electrode). Well 0 [ ] ff ct RM → 0 [ ] ff dl C nF → 0 [ ] ff s Rk → 0 [ ] ff p Hz → 0 [ ] ff z kHz → i ff 1 0.61 2.66 1.27 98.56 32.19 0.02 3 0.44 3.71 1.56 98.08 18.81 0.01 4 0.34 4.76 1.88 97.29 12.28 0.02 5 0.72 2.56 2.01 86.37 21.20 0.02 7 0.35 4.17 1.73 107.95 15.22 0.02 8 0.26 4.78 1.43 129.08 16.05 0.02 Table 3.2 Model parameters estimated at the beginning of the RT simulation (𝑓𝑓 → 0) for AA8 cell line (FO REW model with 10 microelectrodes and one reference electrode). Well 0 [ ] ff ct RM → 0 [ ] ff dl C nF → [ ] s Rk 0 [ ] ff p Hz → 0 [ ] ff z kHz → 1  i ff 1 2.34 4.21 1.71 16.16 15.15 0.95 0.01 3 2.67 6.60 1.44 9.04 14.69 0.93 0.02 4 0.33 4.92 1.56 98.02 14.25 0.99 0.02 5 1.64 3.56 2.09 27.20 14.64 0.96 0.01 7 0.45 5.75 0.71 61.13 26.59 0.96 0.03 8 0.32 5.94 1.16 83.21 15.86 0.97 0.02 98 A similar situation appears with the parameters of the FO REW model of the Table 3.2 if they are compared with the Table 2.12. The values obtained for 𝑓𝑓 → 0 in the RT simulation present less spread than those obtained during the modeling. Fig. 3.4 and Fig. 3.5 illustrate the evolution of ffsim and ffi during the RT simulation of a RT simulation experiment for AA8 cells using the REW and FO REW models, respectively. As can be seen in the figures, the accuracy is considerably high, considering that the values of𝑓 𝐶𝐸 𝑓𝑓→1 and 𝑎𝐶𝐸 𝑓𝑓→1 are not known. Note also that the simulation performed with the REW model seems to be more accurate than the one performed with the FO REW model. The mean relative error (erm) is calculated using the following expression, ( ) ( ) () sim i rm ji ff j ff j eff j − = (3.5) where the sum of j is the sum of the expression for each calculated value of ff. Table 3.3 shows the mean relative errors per well using both models. This table confirms what was previously observed in the Fig. 3.4 and Fig. 3.5, by simulating in RT, more accurate results are obtained, in terms of ff, using the REW (IO) model than the FO REW model. Considering that the FO REW model should fit better to the behavior of the real CE system, the most likely explanation for its lower accuracy is that it is more complex, implying the estimation of more parameters, which can greatly change the other parameters. These parameters that are introduced in the FO REW model are the FO coefficients α1 and α2, which with a small change produce a considerable change in the other parameters, since in terms of bode diagrams (BD) they configure the slope of the magnitude and phase lines, as well as the degrees per pole or zero. The final target of this work is the simulation of a RT experiment to know the cell concentration of a CC assay at each moment. In the following figures and tables, the level of accuracy of the cell concentration will be analyzed compared to the traditional optical counting method. Fig. 3.6 shows the cell concentration estimated in the RT simulation for Nini of 2500, 5000 and 1000 cells. Lines blue and red present the cell concentration obtained by using the FO REW model using ffsim and ffi respectively. Lines yellow and purple present the cell concentration obtained by using the REW model using ffsim and ffi respectively. The last line, in green, present the cell concentration obtained using the traditional optical counting method. As can be seen, as for the ff curves (Fig. 3.4 and Fig. 3.5) the FO REW model performs a little worse than the REW model on the final point (t Table 3.3 Mean relative error (erm) per well for AA8 cell line. erm is shown for REW and FO REW models. Well 1 3 4 5 7 8 [%] REW rm e 11.44 5.76 5.72 11.79 5.11 6.26 [%] FO REW rm e 16.23 6.79 24.73 30.41 9.98 17.65 99 a) b) c) Fig. 3.4 Fill-factor evolution for each well of the experiment using the REW model in the RT simulation. ffsim in blue and ffi in red a) Nini: 2500 cells (well 1 on the left and well 3 on the right). b) Nini : 5000 cells (well 4 on the left and well 5 on the right). c) Nini : 10000 cells (well 7 on the left and well 8 on the right). a) b) c) Fig. 3.5 Fill-factor evolution for each well of the experiment using the FO REW model in the RT simulation ffsim in blue and ffi in red a) Nini: 2500 cells (well 1 on the left and well 3 on the right). b) Nini : 5000 cells (well 4 on the left and well 5 on the right). c) Nini : 10000 cells (well 7 on the left and well 8 on the right). 100 Fig. 3.6 AA8 cell line cell concentration comparison between RT simulation using FO REW model (obtained from ffsim in blue and ffi in red) and REW model (obtained from ffsim in yellow and ffi in purple), and traditional counting method (green ). 101 = 120 h). In addition, the graphs also show that the higher the initial concentration, the more accurate the cell concentration estimation method. As a final point, notice that the RT experiment estimation method calculates the fill factor (and consequently the cell concentration) better at the beginning of the experiment than the estimation of the model parameters carried out in the Chapter 2. Therefore, it is concluded that the estimation of all parameters for each measurement of 𝑓𝐶𝐸 and 𝑎𝐶𝐸 is more accurate than using the parameter values calculated at the end of the experiment. Finally, the mean relative error of the cell concentration (erm) will be calculated using the formulas as follows, . . ( ) ( ) () ( ) ( ) () trad trad sim trad trad trad rm sim jtrad trad i trad trad trad rm i jtrad trad C j C j eCj C j C j eCj − = − =   (3.6) where jtrad is the index of the cell concentration obtained by traditional optical method (defined as Ctrad), Csim is the cell concentration calculated using ffsim, and Ci is the cell concentration estimated using ffi. Since the cell concentration vectors obtained in the RT simulation have many more points than the one obtained by the traditional method (sampling times of 1 h versus 1 day respectively), the error is estimated using the point j closest in time to the jtrad point. Table 3.4 contains the erm calculated from Csim and Ci for the REW and FO REW model. The first point to note is that the average erm is high. However, note that the Csim and Ci curves are closer to or in the error margins of Ctrad. Moreover, the tendency is similar. Therefore, the values are acceptable. On the other hand, although the graphs show that the FO REW model is less accurate at the end than the REW for the RT simulation, the erm values are lower for the FO REW model. The differences between erm.sim and erm.i are almost negligible considering the erm values. Although the erm values of Table 3.4 are better for FO REW model, the graph shows that at the end the REW model gives better results. Finally, note that the FO REW model returns better results in terms of the trend of the curve (first hours of the experiment), as it is similar to that obtained by traditional methods. This may be due to Table 3.4 Cell concentration mean relative error (erm) per Nini for AA8 cell line. erm is shown for REW and FO REW models. The erm obtained from Csim and Ci is calculated for each model. Nini 2500 5000 10000 . [%] REW rm sim e 59.37 37.56 46.40 . [%] REW rm i e 68.47 46.70 49.79 . [%] FO REW rm sim e 33.34 36.79 40.70 . [%] FO REW rm i e 41.04 31.61 31.19 102 the greater complexity of the model, which better describes the real dynamics of the growth of a CC assay. In the conclusions of Chapter 2, it was discussed that the electrical model of the CE block can be fitted using the OBT values of 𝑓𝐶𝐸 and 𝑎𝐶𝐸 at the beginning (𝑓𝑓 → 0) and at the end (𝑓𝑓 → 1) of the experiment. In the current chapter, it has been found that, at least for the best performing cell line (AA8), it is possible to determine the cell concentration in RT during a laboratory experiment. By implementing the algorithm in the OBT prototype (which includes a microcontroller [47]), it is possible to calculate the cell concentration in RT (each time a measurement of𝑓𝐶𝐸 and 𝑎𝐶𝐸 is taken) in an experiment with a CC assay. This is a significant improvement in the field of biology, because the time and economic savings are very large. It is one of the main contributions of this thesis, with which biological and medical research processes can be optimized, and go further with the research or invest the saved resources in another research. 3.1.2.2 N2a results In this section the results of the RT simulation for the N2a line will be presented. The simulation is performed in the same way as for the AA8 cell line (3.1.2.1). Then, the first data to report are the initial parameters derived during the first hours of the RT simulation. Table 3.5 and Table 3.6 shows the initial parameters from simulations using both models, REW and FO REW model. As is the case for the AA8 cell line, the 𝑅𝑐𝑡 𝑓𝑓→0 and 𝑝𝑓𝑓→0 parameters have a higher dispersion using the FO REW model. This is mainly due to the estimation of 𝑝𝑓𝑓→0, which directly influences the calculation of 𝑅𝑐𝑡 𝑓𝑓→0 (equation (2.78)). The values of 𝑝𝑓𝑓→0, in some cases, are much lower (almost an order of magnitude) for the FO REW model than for the REW model. A small fluctuation of α1 can produce these changes, since it changes the slope of the BD of the model. Moreover, taking into account that this cell line has a low 𝛥𝑎𝐶𝐸, and knowing that the value of 𝑝𝑓𝑓→0 is largely determined by 𝑓 𝐶𝐸 𝑓𝑓→0 and 𝑎𝐶𝐸 𝑓𝑓→0 [55], it is expected that more deviations can be seen in the FO REW model as it is more sensible. Table 3.5 Model parameters estimated at the beginning of the RT simulation (𝑓𝑓 → 0) for N2a cell line (REW model). Well 0 [ ] ff ct RM → 0 [ ] ff dl C nF → 0 [ ] ff s Rk → 0 [ ] ff p Hz → 0 [ ] ff z kHz → i ff 2 0.41 3.95 1.68 98.88 16.53 0.004 3 0.38 3.86 1.66 107.94 17.04 0.006 4 0.36 3.69 1.68 121.68 17.67 0.006 6 0.32 4.87 1.43 102.68 15.67 0.004 7 0.32 4.66 1.74 106.60 13.54 0.001 8 0.30 5.49 2.00 95.19 10.00 0.040 103 Table 3.7 displays the erm of the ffsim using as reference ffi. As can be seen, erm in terms of ff is much larger for N2a cell line than for AA8. This was expected, since worse results were obtained in Chapter 2 for the N2a and N2aAPP lines than for AA8, as explained in 2.1.5.5.2. Notice that the erm of the ff does not follow a pattern as a function of Nini. An interesting fact is that the erm of the ff generated from the FO REW model in the simulation is larger than the erm of the ff generated from the REW model. The erm values shown in the table above are obtained from Fig. 3.7 and Fig. 3.8. Figures clearly show very large differences between ffsim and ffi. The hypothesis to explain such high errors is the same as in 2.1.5.5.2, the differences in 𝛥𝑎𝐶𝐸. And the solution, in the same way, is still to tune the BPF to a higher frequency (around 10 kHz). Table 3.6 Model parameters estimated at the beginning of the RT simulation (𝑓𝑓 → 0) for N2a cell line (FO REW). Well 0 [ ] ff ct RM → 0 [ ] ff dl C nF → [ ] s Rk 0 [ ] ff p Hz → 0 [ ] ff z kHz → 1  i ff 2 1.57 6.56 1.58 15.47 10.49 0.94 0.007 3 0.48 4.30 1.88 77.93 13.52 0.99 0.010 4 1.09 6.70 1.28 21.75 12.73 0.93 0.020 6 0.37 4.97 1.81 87.05 12.15 0.99 0.001 7 0.38 5.48 1.45 76.30 13.74 0.98 0.017 8 0.23 5.04 1.45 137.59 15.00 1.01 0.044 Table 3.7 Mean relative error (erm) per well for N2a cell line. erm is shown for REW and FO REW models. Well 2 3 4 6 7 8 [%] REW rm e 19.51 19.00 12.68 7.73 11.87 13.48 [%] FO REW rm e 28.13 52.81 19.29 29.55 33.23 39.82 Table 3.8 Cell concentration mean relative error (erm) per Nini for N2a cell line. erm is shown for REW and FO REW models. The erm obtained from Csim and Ci is calculated for each model. Nini 2500 5000 10000 . [%] REW rm sim e 372.52 154.33 83.80 . [%] REW rm i e 430.92 173.84 91.49 . [%] FO REW rm sim e 205.13 100.82 56.00 . [%] FO REW rm i e 275.54 129.35 55.11 104 While the Fig. 3.9 presents the evolution of Csim and Ci versus Ctrad for the N2a cell line, the Table 3.8 shows the errors of these concentrations. Of course, the Ctrad time vector had to be advanced 48 h to match the dynamics, as discussed in 2.1.5.5.1. The figure shows that, in terms of the final value of cell concentration (t = 120 h), Csim is always more distant than Ci from Ctrad. This does not necessarily mean that it is a worse fit, since the table shows slightly worse erm values for Ci than for Csim for any Nini. On the other hand, note that Ci in the figure fits much better to Ctrad if the first hours of the simulation are excluded. Finally, although the cell concentration obtained from the FO REW model achieves significantly lower erm, notice that it is closer to the real dynamics at the beginning of the experiment but from the middle to the end of the simulation, the FO REW model does not produce as good results as the REW model. 105 a) b) c) Fig. 3.7 Fill-factor evolution for each well of the experiment using the REW model in the RT simulation. ffsim in blue and ffi in red a) Nini: 2500 cells (well 2 on the left and well 6 on the right). b) Nini : 5000 cells (well 3 on the left and well 7 on the right). c) Nini : 10000 cells (well 4 on the left and well 8 on the right). a) b) c) Fig. 3.8 Fill-factor evolution for each well of the experiment using the FO REW model in the RT simulation ffsim in blue and ffi in red a) Nini: 2500 cells (well 1 on the left and well 3 on the right). b) Nini : 5000 cells (well 4 on the left and well 5 on the right). c) Nini : 10000 cells (well 7 on the left and well 8 on the right). 112 second one providing much better results (at the accuracy level of FO REW model). Another point to note is that, in general, better results are obtained at higher Nini. Finally, it should be pointed out that the first section of the simulation is not as important in terms of predicting ff and cell concentration. Therefore, these data are quite acceptable, and provide a useful starting point for an improvement in the models and the method of simulation and parameter computation. Separate analysis of the two zones of the experiment has provided interesting results. The cell concentration of a CC assay can be estimated in RT by connecting the CC assay as a load to an OBT. The errors are still large, but with some improvements of the algorithm, and of the OBT measurement prototype, the error should be greatly reduced. 3.2 Electrical simulation of the Bio-Oscillator In the previous sections and chapters, the cell concentration in a CC assay has been obtained from the oscillation data of the electronic system where the CE block is connected as a load to the OBT circuit. Apparently, the results are acceptable as the obtained cell concentration has been compared with the measured cell concentration obtained by traditional optical counting method. Oscillation amplitude and frequency are the key parameters to achieve the main goal: provide the cell concentration of a CC assay experiment in RT. On the other hand, the oscillator must be designed to produce self-sustained oscillations. This point requires to perform electrical simulations (ES) of the oscillator, which should consider similar load values to the CE block. In this way it is possible to know the response of the system to different values of ff and other parameters of the electrical model. Prior to performing the experiments to obtain the experimental amplitude and frequency values, electrical simulations of the oscillator were carried out considering values of the model parameters obtained theoretically. In future work, the oscillator prototype will be modified, so it is crucial to have very accurate values of the model parameters, so that the prototype and the results obtained are as optimal and accurate as possible. Therefore, because of the relevance of the electrical simulations, it should be ensured that the values of the model parameters calculated in the previous sections can be used in future electrical simulations. In the current subsection ES of the oscillator will be performed, using the values of the model parameters previously obtained. The simulations will be performed using the electrical simulator Ngspice [61], but with the difference that all the electrical simulations are launched from Matlab. As a comparison metric, the amplitude and frequency data will be obtained using the method of minimization of a CF, and the error with respect to the amplitudes and frequencies obtained in a real experiment will be calculated. The same error will be estimated for the results of the electrical simulations and the error of both methods will be compared. As in the previous subsections, tests will be performed using data from cell lines AA8, N2a and N2aAPP. 113 3.2.1 Electrical simulator and its integration in Matlab The electrical simulator used is Ngspice, which is an open-source Spice simulator for electric and electronic circuits. All its documentation can be found at [61]. Ngspice is used in this work due to the fact that it can be launched from Matlab using commands through the system function. To simulate the oscillator circuit in order to reproduce each point of the experiments performed, it is necessary to perform as many simulations as amplitude and frequency measurements have been taken in each of the experiments with the different cell lines. Considering that the time between each sample taken is 1 h, and that the experiments have a duration of between 6 and 9 days, at least 170 simulations must be performed for each cell line. This is unmanageable, since editing each of the components of the electrical model of the CE block one by one would be too laborious. All this without considering that exporting the results to Matlab for processing becomes really impractical. For these reasons, Matlab and Ngspice are connected to each other to edit, simulate and process the data automatically. Section 3.1.1 presents the simulation routine of the system under test using mathematical methods of a CF minimization. These simulations emulate a RT experiment in which the future values of 𝑓𝐶𝐸 and 𝑎𝐶𝐸 are not known, recalculating certain model parameters for each measurement of f and a. In this way, acceptable results were achieved in predicting the behavior of the cell concentration of a CC assay. The initial parameters calculated in mentioned section, together with the parameters that are reevaluated for each measurement of 𝑓𝐶𝐸 and 𝑎𝐶𝐸, will be introduced in the electrical simulator to reproduce the frequency and amplitude values obtained experimentally. The whole circuit (oscillator plus the electrical model of the CE block has been implemented in Spice language with the same topology and characteristics as described in section 2.1.4. Notice that the file which contains the electrical model of the CE block can be easily replaced to switch between the different electrical models (SEW, REW and FO REW). The details, step by step, of the electrical simulation using Matlab in combination with Ngspice are described below: 1. Edit the model netlist. The oscillator circuit is not modified from one measurement to another. Therefore, the part of the circuit to be edited is the CE block model. This block, composed of passive components, is defined by the model used, which can be SEW, REW and FO REW. There is a text file in Spice language defining the structure of the CE model to be used in the simulations. For each actual measurement of 𝑓𝐶𝐸 and 𝑎𝐶𝐸, the CF minimization method obtained the values of the model parameters. Matlab edits the Spice file that contain the model to include the values of the passive components of the CE model. Matlab does it automatically, using string replacement. In this way it is not necessary to edit the values one by one before each electrical simulation. 114 2. Run the simulation. As previously mentioned, an electrical simulation with Ngspice can be executed from Matlab using commands. The simulation is run in a single command line, choosing certain options, the Spice file to be simulated and the files where the results are stored. Using commands, it is not possible to choose the analysis time to be performed (DC, transient, AC, etc.), neither the variables to be displayed and saved at the end of the simulation. All this is included in the Spice file to be simulated. For this reason, a script has been designed to which the Spice file to be simulated, the analysis to be performed, the necessary parameters for the analysis and the variables to be saved at the end of the simulation are provided. This Matlab script or function works as a mask to automate the process of editing the type of analysis and the results to be shown, since for a user not used to the use of Spice, it can be useful to launch in a simple way from Matlab some electrical simulations in Spice (for example, for educational purposes). As in step 1, the script inserts the necessary code into the Spice file to perform the required analysis. 3. Results generation. After finishing the simulation, the script described in step 2 extracts the variables to be saved from a text file generated by Ngspice. This file contains the information of the variables that have been requested to be saved. These variables are stored in a Matlab struct variable that contains their names and values, as well as other variables that may be necessary for their interpretation (e.g., time). Using the mixed electrical simulation technique which involves Matlab and Ngspice, it is possible to automate the process of CE block configuration, electrical simulation and data analysis, which is much more efficient than performing the simulations one by one. All the scripts to perform the simulations using Matlab are in Appendix B. The function simOBTngspice is the main function, managing the whole simulation. Function editBioZngspice edits automatically the netlist of the CE model used. ngspiceSim is the function that manages the simulation of a spice circuit, launching Ngspice from Matlab and processing the results. The function getfA obtains the frequency and amplitude of a sine wave without using the Fast Fourier Transform. 3.2.2 Electrical simulation using REW electric model Due to the fact that the SEW model is less accurate, responds worse to possible perturbations (an example is the values of 𝑓𝐶𝐸 and 𝑎𝐶𝐸 at the beginning of an experiment) and the parameter values obtained are not close to those theoretically expected (because it is not the real model), it is decided to implement the REW model for the electrical simulations. These electrical simulations are performed as explained in the previous subsection and using the data from the three cell lines. As an example of the result of the electrical simulation, Fig. 3.13 shows the 𝑉𝐶𝐸voltage waveform (at the bottom of the 115 figure) obtained in that simulation. A waveform is obtained for each time of the simulation (for each t(j)), from which values of 𝑓𝐶𝐸 and 𝑎𝐶𝐸 are extracted for each t(j) (top of the figure). The value of Rsi obtained at the beginning of the RT simulation (AA8 in Table 3.1, N2a in Table 3.5 and N2aAPP in Table 3.9), and the values of ff, Rct, Cdl, ΔRs and Rgap obtained for each measure of f and a are used. Since there is a value of each parameter for each measure of 𝑓𝐶𝐸 and 𝑎𝐶𝐸, they are not shown because of the space requirement in this document. In addition, using the CF minimization method, the values of 𝑓𝐶𝐸 and 𝑎𝐶𝐸 will be estimated from these same parameters. As a comparison metric, the erm in terms of frequency and amplitude is used in relation to the real values measured experimentally for each cell line. erm is calculated as follows, 100 100 f CEmeas CEestim rm CEmeas a CEmeas CEestim rm CEmeas ff ef aa ea − = − = (3.7) where 𝑒𝑟𝑚 𝑓 and 𝑒𝑟𝑚 𝑎 are the mean relative error in terms of frequency and amplitude respectively, 𝑓𝐶𝐸𝑚𝑒𝑎𝑠 and 𝑎𝐶𝐸𝑚𝑒𝑎𝑠 are the values of 𝑓𝐶𝐸 and 𝑎𝐶𝐸, respectively, measured in a real experiment, and 𝑓𝐶𝐸𝑒𝑠𝑡𝑖𝑚 and 𝑎𝐶𝐸𝑒𝑠𝑡𝑖𝑚 are the values of 𝑓𝐶𝐸 and 𝑎𝐶𝐸, respectively, obtained by one of the both methods (CF minimization and electrical simulation). Amplitude Frequency Fig. 3.13 Voltage 𝑉𝐶𝐸 waveform (bottom graph) obtained from a Ngspice electrical simulation of the oscillator circuit (including the CE block).From the wave,𝑓𝐶𝐸 and 𝑎𝐶𝐸 are extracted for each time t(j) (top graphs). 116 In the following subsections, the frequency and amplitude erm data for the three cell lines obtained using the REW model will are presented. 3.2.2.1 AA8 In this subsection the erm values, in terms of frequency and amplitude, obtained with CF minimization and by electrically simulating the oscillator using REW model and AA8 cell line data are presented. In Table 3.15 and Table 3.16 erm values are shown. 𝑒𝑟𝑚.𝑓 𝐶𝐹 and 𝑒𝑟𝑚.𝑎 𝐶𝐹 are the mean relative errors of the frequency and amplitude estimated by using the CF minimization method. On the other hand, 𝑒𝑟𝑚.𝑓 𝐸𝑆 and 𝑒𝑟𝑚.𝑎 𝐸𝑆 are the mean relative errors of the frequency and amplitude obtained from the electrical simulation (ES). The frequency error is smaller than the amplitude error, which may be due to the smaller amplitude ranges. On the other hand, as expected, the error using the CF minimization method is lower than that obtained for the ESs. Larger errors were expected for ES because the model parameter values are calculated for CF minimization, not for ES. Therefore, it is reasonable to expect bigger errors due to tolerances and computational methods. For the AA8 cell line, it is concluded that the errors found during the ES are acceptable. In the future, simulations will be performed based on this data to redesign or improve the oscillator prototype, since the errors are acceptable and do not imply a problem. Perhaps, if the circuit were to be integrated, more accurate data would be needed, since smaller voltage ranges would be used. Table 3.15 Mean relative error (erm) per well for AA8 cell line in terms of frequency. Well 1 3 4 5 7 8 . [%] CF rm f e 0.34 0.31 0.27 0.25 0.28 0.16 . [%] ES rm f e 1.34 0.55 1.84 1.77 1.30 1.26 Table 3.16 Mean relative error (erm) per well for AA8 cell line in terms of amplitude. Well 1 3 4 5 7 8 . [%] CF rm a e 1.14 1.25 1.24 0.82 1.01 0.84 . [%] ES rm a e 6.36 9.14 8.62 6.46 8.55 7.85 117 3.2.2.2 N2a As in the previous subsection, the erm values, in terms of frequency and amplitude, obtained with CF minimization and by electrically simulating the oscillator using REW model and N2a cell line data are presented. Table 3.17 and Table 3.18 show the values of the erm for the N2a cell line. The conclusions reached for the error data of the AA8 cell line are applicable to the N2a cell line. In addition, it should be added that the errors of the N2a line are slightly higher than the AA8. This fact may be due, as explained in previous sections, to the difference in the ranges of f and a between the different cell lines. However, the errors are also acceptable, which means that these data can be used in future ESs. 3.2.2.3 N2aAPP The erm values, in terms of frequency and amplitude, obtained from CF minimization and electrical simulation using REW model and N2aAPP cell line data are presented in the current subsection. Table 3.19 and Table 3.20 show the frequency and amplitude erm values respectively. All comments and conclusions given for the two previous cell lines are applicable to the N2aAPP cell line. The error of both parameters is similar to the error values of the previous cell lines. All three cell lines have acceptable error values. The Table 3.17 Mean relative error (erm) per well for N2a cell line in terms of frequency. Well 2 3 4 6 7 8 [%] CF rm e 0.51 0.36 0.30 0.37 0.28 0.33 [%] ES rm e 2.21 1.97 1.65 2.08 2.01 2.47 Table 3.18 Mean relative error (erm) per well for N2a cell line in terms of amplitude. Well 2 3 4 6 7 8 [%] CF rm e 1.94 1.59 1.48 1.55 1.27 1.39 [%] ES rm e 10.58 9.06 9.25 9.70 8.84 10.11 Table 3.19 Mean relative error (erm) per well for N2aAPP cell line in terms of frequency. Well 2 3 4 6 7 8 [%] CF rm e 0.39 0.36 0.37 0.25 0.24 0.90 [%] ES rm e 2.20 1.69 1.58 1.67 1.76 1.09 118 values of the REW model parameters obtained are valid for performing ESs, which allows testing of future improvements of the prototype. 3.3 Conclusions The main target of Chapter 3 is to simulate in RT a CC assays growth real experiment connected to an OBT as a load and make predictions. The simulation has the goal to estimate the cell concentration in the wells in RT, i.e., when a measurement of 𝑓𝐶𝐸 and 𝑎𝐶𝐸 is performed, calculating the cell concentration at that time. In order to perform this task, the parameters of the electrical model used to predict the electrical behavior of the CE block must be estimated. Keeping in mind that future values are not known, so parameters such as Rgap, which has more influence at the end of the simulation (𝑓𝑓 → 1), must be calculated for each measurement taken. The method of minimization of a CF has been used to obtain the optimal values of the model parameters that meet the BSC for each value of 𝑓𝐶𝐸 and 𝑎𝐶𝐸. Due to the complexity of the used models (REW and FO REW) it is hard to obtain good results, since there is a considerable number of parameters to calculate at each moment, and small variations in one of them can make the value of the parameter ff differs from its real value. In addition, the 𝑓𝐶𝐸 and 𝑎𝐶𝐸 data for some cell lines are worse than the data for others, as the prototype measurement is still in the experimental phase. Even with these difficulties, it has been possible to successfully estimate the RT cell concentration present in a CC assay, although with a certain margin of error. Considering the variability of cell concentration, and the studies performed up to now, it can be considered that the results are good. It is a very useful method, but should be improved in future works. Some conclusions have been derived from these simulations, which are discussed below: • Cell concentration of a CC assay can be measured using an oscillator as a measuring circuit to which the CC is connected as a load. The oscillator circuit (OBT) must be configured to oscillate at the frequencies at which the most robust measurements are obtained taking into account the frequency response (BD) of the CC assay (CE block). • The optimum oscillation frequency depends on the cell line. In this work, all measurements have been performed with the cut-off frequency of the BPF of the OBT set at 1 kHz, i.e., the circuit oscillates at frequencies somewhat lower than the filter cut-off frequency. The cell lines used as subject of study in this Table 3.20 Mean relative error (erm) per well for N2aAPP cell line in terms of amplitude. Well 2 3 4 6 7 8 [%] CF rm e 1.43 1.53 1.68 0.93 0.90 3.33 [%] ES rm e 9.58 8.60 9.06 9.44 10.25 8.21 119 work (AA8, N2a and N2aAPP) have provided relatively good results in terms of oscillation parameters. However, after observing the BD of other different cell lines, in future works, the cut-off frequency of the BPF should be increased to obtain larger ranges of frequency and amplitude (𝛥𝑓𝐶𝐸 and 𝛥𝑎𝐶𝐸). • The erm obtained, in terms of cell concentration, depends on the cell line, because the lines with a lower Acell (N2a and N2aAPP) reach much larger error values than the cell line with a higher Acell (AA8). • More generally, it has been demonstrated that the cell concentration can be obtained in RT during the performance of a cell growth experiment. The errors found, while significant, can be reduced by improving the measurement circuit and the algorithm for calculating the cell concentration. The other purpose of Chapter 3, derived from the main goal, is to obtain the values of the parameters of the CE model of the CE block. Although it is considered secondary in this work, it is crucial for the development of the measurement prototype, since it allows to perform ES with real data of the CE model. The ES are very important to improve and test the OBT before implementing improvements to the current prototype or designing a new one. For this reason, ES have been performed with the CE model data using the Ngspice electrical simulator. Routines have been created in Matlab to automatically edit the Spice files of the CE model and to launch simulations with Ngspice in an automated and much easier way than using the Ngspice software's own command window. The results of the electrical simulations are very promising, as the error margins are acceptable, so the CE model data obtained can be used as a starting point to improve the OBT circuit and obtain more robust and accurate measurements of the RT cell concentration of a CC assay. 120 121 Chapter 4 Conclusions and future works The goal of this thesis is to find a method to estimate in real-time (RT) the cell concentration in a cell culture (CC) assay real experiment using bioimpedance (BI) measurements. For this purpose, an electronic oscillator circuit (OBT) has been used which ensures the measurement of BI in a non-invasive, fast and safe way for the cells. Some conclusions have been reached after finishing this thesis: • The BI of a biological sample (BS) is usually obtained by injecting a signal of a specific frequency and amplitude through the BS. This is a limitation, because although it is possible to calculate the impedance of the BS at a particular frequency or frequencies, it is very complex to obtain an electrical model of these averages. In this thesis an indirect BI measurement technique has been used. The CC is connected to an OBT circuit, which oscillates by itself at an amplitude and frequencies determined by all the elements that compose the circuit, including the CC. The oscillation frequency and amplitude depend, therefore, on the BI of the CC. From the two oscillation parameters, it is possible to obtain the electrical model of the cell-electrode (CE) block. • Although the models tested in this work consider the contribution of the cells to the electrical behavior of the CE block to be negligible (not including capacitance or intracellular resistance), these models are reasonably accurate at predicting the electrical response of the block by selecting the adequate frequency. This is because at some frequencies the cells impedance achieves very high values (short circuit). The fill-factor (ff) parameter can measure, with high precision, the fill percentage of the bottom of the well, which means that the cell concentration can be obtained in an approximate range. This means that it is possible to perform a RT monitoring of the CC assay growth from ff = 0 to ff = 1, which is a key contribution of this work (only applicable to a monolayer CC assay). • The electrode model used has been widely employed in the BI field, and in this work, it has provided good results. Although some parameters differ a little from the theoretical calculated values, this is not a problem, because there are many variables during the theoretical calculation whose values are still under study. 128 Sensors, vol. 14, no. 6, pp. 10895–10928, Jun. 2014. [15] Y.-Y. Lu, J.-J. Huang, and K.-S. Cheng, “The design of electrode-array for monitoring the cellular bioimpedance,” in 2009 {IEEE} Symposium on Industrial Electronics & Applications, Oct. 2009. [16] K. Lei, “Review on Impedance Detection of Cellular Responses in Micro/Nano Environment,” Micromachines, vol. 5, no. 1, pp. 1–12, 2014. [17] D. A. Borkholder, “Cell based biosensors using microelectrode,” Stanford University, 1998. [18] R. Pradhan, M. Mandal, A. Mitra, and S. Das, “Monitoring cellular activities of cancer cells using impedance sensing devices,” Sensors Actuators B Chem., vol. 193, pp. 478–483, Mar. 2014. [19] P. Daza, A. Olmo, D. Cañete, and A. Yúfera, “Monitoring living cell assays with bio-impedance sensors,” Sensors Actuators B Chem., vol. 176, pp. 605–610, 2013. [20] M. Abdolahad, H. Shashaani, M. Janmaleki, and S. Mohajerzadeh, “Silicon nanograss based impedance biosensor for label free detection of rare metastatic cells among primary cancerous colon cells, suitable for more accurate cancer staging,” Biosens. Bioelectron., vol. 59, pp. 151–159, Sep. 2014. [21] C. F. Lourenço, A. Ledo, J. Laranjinha, G. A. Gerhardt, and R. M. Barbosa, “Microelectrode array biosensor for high-resolution measurements of extracellular glucose in the brain,” Sensors Actuators B Chem., vol. 237, pp. 298– 307, 2016. [22] A. Dibao-Dina, J. Follet, M. Ibrahim, A. Vlandas, and V. Senez, “Electrical impedance sensor for quantitative monitoring of infection processes on {HCT}-8 cells by the waterborne parasite Cryptosporidium,” Biosens. Bioelectron., vol. 66, pp. 69–76, 2015. [23] S. Reitinger, J. Wissenwasser, W. Kapferer, R. Heer, and G. Lepperdinger, “Electric impedance sensing in cell-substrates for rapid and selective multipotential differentiation capacity monitoring of human mesenchymal stem cells,” Biosens. Bioelectron., vol. 34, no. 1, pp. 63–69, 2012. [24] S. Grimnes and O. G. Martinsen, Bioimpedance and Bioelectricity Basics, 3rd Ed. Oslo: Elsevier, 2015. [25] X. Huang, D. Nguyen, D. W. Greve, and M. M. Domach, “Simulation of Microelectrode Impedance Changes Due to Cell Growth,” {IEEE} Sensors J., vol. 4, no. 5, pp. 576–583, Oct. 2004. [26] A. Olmo and A. Yúfera, “Computer simulation of microelectrode based bioimpedance measurements with COMSOL,” in Third International Conference on Biomedical Electronics and Devices, 2010, pp. 178–182. [27] X. Zhang, W. Wang, and S. Jang, “The Modeling, Design, Fabrication, and Application of Biosensor Based on Electric Cell-Substrate Impedance Sensing (ECIS) Technique in Environmental Monitoring,” Biosens. Environ. Monit., Dec. 2018. [28] I. Giaever and C. R. Keese, “Micromotion of mammalian cells measured electrically,” Proc. Natl. Acad. Sci. U. S. A., vol. 88, no. 17, pp. 7896–7900, 1991. [29] G. T. A. Kovacs, “Microelectrode models for neural interfaces,” in Enabling 129 technologies for cultured neural networks, D. A. Stenger and T. M. McKenna, Eds. New York: Academic Press, 1994, pp. 121–165. [30] J. A. Serrano, P. Pérez, A. Maldonado, M. Martín, A. Olmo, P. Daza, G. Huertas, and A. Yúfera, “Practical characterization of cell-electrode electrical models in bio-impedance assays,” in BIODEVICES 2018 - 11th International Conference on Biomedical Electronics and Devices, Proceedings; Part of 11th International Joint Conference on Biomedical Engineering Systems and Technologies, BIOSTEC 2018, 2018, vol. 1, pp. 100–108. [31] A. Yúfera, P. Daza, and D. Cañete, “Using microelectrode models for real time cell-culture monitoring,” in Proceedings of the Annual International Conference of the IEEE Engineering in Medicine and Biology Society, EMBS, 2011, vol. 2011, pp. 3983– 3986. [32] F. Liu, S. M. Arifuzzaman, A. N. Nordin, D. Spray, and I. Voiculescu, “Characterization of endothelial cells using electrochemical impedance spectroscopy,” IEEE Asia-Pacific Conf. Circuits Syst. Proceedings, APCCAS, pp. 252– 255, 2010. [33] N. Goda, N. Kataoka, J. Shimizu, S. Mohri, Y. Yamamoto, and F. Kajiya, “Evaluation of micromotion of vascular endothelial cells with estrogen in electrical cell-substrate impedance sensing (ECIS) method using a mathematical model - Cell-to-cell distance and cell-to-substrate distance,” 2005 3rd IEEE/EMBS Spec. Top. Conf. Microtechnology Med. Biol., vol. 2005, pp. 211–214, 2005. [34] J. J. Cabrera, “Electrical Characterization of Vegetable Tissues Based on Fractional- Order Models,” 2020. [35] A. Maldonado, P. Perez, G. Huertas, A. Yufera, A. Rueda, and J. L. Huertas, “Monitoring tissue evolution on electrodes with bio-impedance test,” 2016 Conf. Des. Circuits Integr. Syst. DCIS 2016 - Proc., pp. 213–218, Feb. 2017. [36] J. O. Bockris and A. K. N. Reddy, Modern Electrochemistry: Volume 1. Boston, MA: Springer US, 1970. [37] D. A. Robinson, “The Electrical Properties of Metal Microelectrodes,” Proc. IEEE, vol. 56, no. 6, pp. 1065–1071, 1968. [38] H. P. Schwan, “Determination of biological impedance,” in Physical in Biological Research, W. L. Nastuk, Ed. New York: Academic Press, 1963. [39] S. R. W, B. J. G, and S. H. P, “Nonlinear AC and DC polarization of platinum electrodes.,” IEEE Trans Biomed Eng, vol. 27, no. 3, pp. 166–171, 1980. [40] H. P. Schwan, “Linear and nonlinear electrode polarization and biological materials,” Ann. Biomed. Eng. 1992 203, vol. 20, no. 3, pp. 269–288, May 1992. [41] J. J. Ackmann, “Complex bioelectric impedance measurement system for the frequency range from 5 Hz to 1 MHz,” Ann. Biomed. Eng. 1993 212, vol. 21, no. 2, pp. 135–146, Mar. 1993. [42] P. R and W. JG, “Bioelectric impedance measurements using synchronous sampling,” IEEE Trans. Biomed. Eng., vol. 40, no. 8, pp. 824–829, 1993. [43] G. Huertas, A. Maldonado, and A. Yúfera, “Electrical Cell-substrate Impedance Spectroscopy (ECIS) measurements based on Oscillation-Based Test techniques,” Lect. Notes Impedance Spectrosc., pp. 109–118, Jan. 2015. 130 [44] R. Bragós, E. Sarró, A. Fontova, A. Soley, J. Cairó, A. Bayés-Genís, and J. Rosell, “Four versus two-electrode measurement strategies for cell growing and differentiation monitoring using electrical impedance spectroscopy,” Annu. Int. Conf. IEEE Eng. Med. Biol. - Proc., pp. 2106–2109, 2006. [45] R. D. Beach, R. W. Conlan, M. C. Godwin, and F. Moussy, “Towards a miniature implantable in vivo telemetry monitoring system dynamically configurable as a potentiostat or galvanostat for two- and three-electrode biosensors,” IEEE Trans. Instrum. Meas., vol. 54, no. 1, pp. 61–72, Jan. 2005. [46] “Applied Biophysics.” http://www.biophysics.com/ (accessed Nov. 05, 2020). [47] P. Pérez, “Circuit design for biomedical laboratories based on impedance measurement,” Universidad de Sevilla, 2019. [48] R. N. Bracewell, The Fourier Transform and Its Applications (3rd ed.). Boston: McGraw-Hill, 2000. [49] J. Wegener, C. R. Keese, and I. Giaever, “Electric cell-substrate impedance sensing (ECIS) as a noninvasive means to monitor the kinetics of cell spreading to artificial surfaces,” Exp. Cell Res., vol. 259, no. 1, pp. 158–166, Aug. 2000. [50] O. Stern, “Zur Theorie der Elektrolytischen Doppelschicht,” Zeitschrift fur Elektrochemie, vol. 30, no. 21‐22, pp. 508–516, Nov. 1924. [51] H. Helmholtz, “Studien über electrische Grenzschichten,” Ann. Phys., vol. 243, no. 7, pp. 337–382, Jan. 1879. [52] M. Gouy, “Sur la constitution de la charge électrique à la surface d’un électrolyte,” J. Phys. Théorique Appliquée, vol. 9, no. 1, pp. 457–468, 1910. [53] D. L. Chapman, “ LI. A contribution to the theory of electrocapillarity ,” London, Edinburgh, Dublin Philos. Mag. J. Sci., vol. 25, no. 148, pp. 475–481, Apr. 1913. [54] H. J. Lewerenz, “On the Structure of the Helmholtz Layer and its Implications on Electrode Kinetics,” ECS Trans., vol. 50, no. 52, pp. 3–20, Apr. 2013. [55] J. Serrano, G. Huertas, A. Maldonado-Jacobi, A. Olmo, P. Pérez, M. Martín, P. Daza, and A. Yúfera, “An Empirical-Mathematical Approach for Calibration and Fitting Cell-Electrode Electrical Models in Bioimpedance Tests,” Sensors, vol. 18, no. 7, p. 2354, Jul. 2018. [56] P. Pérez, A. Maldonado-Jacobi, A. J. López, C. Martínez, A. Olmo, G. Huertas, and A. Yúfera, “Remote Sensing of Cell-Culture Assays,” in New Insights into Cell Culture Technology, InTech, 2017, pp. 135–155. [57] P. Pérez, G. Huertas, A. Maldonado-Jacobi, M. Martín, J. A. Serrano, A. Olmo, P. Daza, and A. Yufera, “Sensing Cell-Culture Assays with Low-Cost Circuitry,” Sci. Rep., vol. 8, no. 1, p. 11, 2018. [58] P. Pérez, G. Huertas, A. Olmo, A. Maldonado-Jacobi, J. Serrano, M. Martín, P. Daza, and A. Yúfera, “Remote Cell Growth Sensing Using Self-Sustained Bio- Oscillations,” Sensors, vol. 18, no. 8, p. 2550, Aug. 2018. [59] G. Huertas, A. Maldonado, A. Yufera, A. Rueda, and J. L. Huertas, “The Bio- Oscillator: A Circuit for Cell-Culture Assays,” IEEE Trans. Circuits Syst. II Express Briefs, vol. 62, no. 2, pp. 164–168, Feb. 2015. [60] “Matlab.” https://es.mathworks.com (accessed Nov. 05, 2020). 131 [61] “Ngspice.” http://ngspice.sourceforge.net (accessed Nov. 05, 2020). [62] “Texas Instruments.” https://www.ti.com/ (accessed Nov. 16, 2020). [63] Y. Yuste, J. A. Serrano, A. Olmo, A. Maldonado-Jacobi, P. Pérez, G. Huertas, S. Pereira, F. de la Portilla, and A. Yúfera, “Monitoring muscle stem cell cultures with impedance spectroscopy,” in BIODEVICES 2018 - 11th International Conference on Biomedical Electronics and Devices, 2018, vol. 1. [64] C. a Monje, Y. Q. Chen, B. M. Vinagre, D. Xue, and V. Feliu, Fractional-order Systems and Controls. Fundamentals and Applications. 2010. [65] “fminsearch.” https://es.mathworks.com/help/matlab/ref/fminsearch.html (accessed Mar. 12, 2021). [66] J. A. Serrano, P. Perez, G. Huertas, and A. Yufera, “Alternative General Fitting Methods for Real-Time Cell-Count Experimental Data Processing,” IEEE Sens. J., vol. 20, no. 24, pp. 15177–15184, Dec. 2020. [67] J. A. Serrano, A. J. Gines, and E. Peralias, “Fast Simulation of Non-Linear Circuits using Semi-Analytical Solutions Based on the Matrix Exponential,” in IEEE International Symposium on Circuits and Systems (ISCAS), 2020. [68] P. Pérez, J. A. Serrano, M. E. Martín, P. Daza, G. Huertas, and A. Yúfera, “A computer-aided design tool for biomedical OBT sensor tuning in cell-culture assays,” Comput. Methods Programs Biomed., vol. 200, p. 105840, Mar. 2021. 132 133 Appendix A: Scientific publications The present thesis has contributed to the publication of some papers in important journals (Table A.1) and the attendance and publication in some international conferences (Table A.2). Table A.1 Journal publications Name Journal Year Quartile Reference Sensing cell-culture assays with low-cost circuitry Scientific reports 2018 Q1 [57] An empirical-mathematical approach for calibration and fitting cell-electrode electrical models in bioimpedance tests Sensors MDPI 2018 Q1 [55] Remote cell growth sensing using self-sustained biooscillations Sensors MDPI 2018 Q1 [58] 3D-printed sensors and actuators in cell culture and tissue engineering: framework and research challenges Sensors MDPI 2020 Q1 [8] Alternative General Fitting Methods for Real-Time Cell- Count Experimental Data Processing IEEE Sensors Journal 2020 Q2 [66] A computer-aided design tool for biomedical OBT sensor tuning in cell-culture assays Computer Methods and Programs in Biomedicine 2021 Q1 [68] Table A.2 Conference publications Name Conference Year Reference Monitoring muscle stem cell cultures with impedance spectroscopy 11th BIODEVICES 2018 [63] Practical characterization of cell-electrode electrical models in bio-impedance assays 11th BIODEVICES 2018 [30] Designing bioimpedance based sensors for cell cultures test 25th DCIS 2020 [6] 134 135 Appendix B: Ngspice toolbox code Main Function function [f,A,R1,R2,C1,C2] = simOBTngspice(R,C,Rgap,Rsi,iRs,ff,n,elec,ref) R1 = R/(1-ff); R2 = R/ff; C1 = C*(1-ff); C2 = C*ff; Rs = Rsi+iRs*(ff^n); if elec==0 ref = []; else ref.Rs = (Rsi+iRs*(ff^n))*sqrt(ref.ker); ref.R1 = ref.Rr/(1-ff); ref.R2 = ref.Rr/ff; ref.C1 = ref.Cr*(1-ff); ref.C2 = ref.Cr*ff; end % Edit circuit out = editBioZngspice (R1,R2,C1,C2,Rgap,Rs,elec,ref); % NGSpice Simulation netlistfilename = 'ngspice_OBT.cir'; simType = {'tran'}; tstep = 1e-005; tfinal = 0.001*32; tini = 0.001*25; maxstep = 1e-005; simParams = {[tstep tfinal tini maxstep]}; outVars = {'0','vvcmp','vvbpf','vvcs','vvpa'}; data_ngs = ngspiceSim(netlistfilename,simType,simParams,outVars); % Estimate amplitude and frequency t = data_ngs.time; 136 L = length(t); y = data_ngs.outVars(3).value; Lint = 2^nextpow2(L); tint = linspace(t(1),t(end),Lint); ts = tint(2)-tint(1); yint = interp1(t,y,tint); ymean = mean(yint); out = getFA(yint-ymean,tint,Lint); f = out.f; A = out.A; end Edit Bioimpedance to Ngspice function out = editBioZngspice (R1,R2,C1,C2,Rgap,Rs,elec,ref) fclose ('all'); % Parameters flagR1 = 'R1val'; flagR2 = 'R2val'; flagC1 = 'C1val'; flagC2 = 'C2val'; flagRgap = 'Rgapval'; flagRs = 'Rsval'; if elec==1 flagR1_r = 'R1val_r'; flagR2_r = 'R2val_r'; flagC1_r = 'C1val_r'; flagC2_r = 'C2val_r'; flagRgap_r = 'Rgapval_r'; flagRs_r = 'Rsval_r'; end % Copy base code of Bioimpedance circuit if elec==0 circuito = ‘ngspice_BioZ_sub.mod'; circuito_base = 'ngspice_BioZ_sub_base.mod'; else 137 circuito = 'ngspice_BioZ11e_sub.mod'; circuito_base = 'ngspice_BioZ11e_sub_base.mod'; end fid = fopen(circuito_base); fulltext = fread(fid,'*char')'; close = fclose(fid); if close==-1 out = close; disp('Error al cerrar durante la lectura del circuito Spice'); return; else out = close; end % Write parameters values fulltext = strrep(fulltext,flagR1_r,num2str(ref.R1)); fulltext = strrep(fulltext,flagR2_r,num2str(ref.R2)); fulltext = strrep(fulltext,flagC1_r,num2str(ref.C1)); fulltext = strrep(fulltext,flagC2_r,num2str(ref.C2)); fulltext = strrep(fulltext,flagRgap_r,num2str(ref.Rgap)); fulltext = strrep(fulltext,flagRs_r,num2str(ref.Rs)); fulltext = strrep(fulltext,flagR1,num2str(R1)); fulltext = strrep(fulltext,flagR2,num2str(R2)); fulltext = strrep(fulltext,flagC1,num2str(C1)); fulltext = strrep(fulltext,flagC2,num2str(C2)); fulltext = strrep(fulltext,flagRgap,num2str(Rgap)); fulltext = strrep(fulltext,flagRs,num2str(Rs)); % Save Spice files fid = fopen(circuito,'w'); fprintf(fid,'%s',fulltext); close = fclose(fid); if close==-1 out = close; disp('Error al cerrar durante la escritura del circuito Spice'); return; else out = close; 144 data(i).plotname = char(regexp(rawfile(simPlotNameIni(i):simFlagsIni(i)), … '^Plotname:\s*(\w.*?\w)\s*$', 'tokens', 'once', 'lineanchors')); data(i).flags = char(regexp(rawfile(simFlagsIni(i):simNoVariablesIni(i)), … '^Flags:\s*(\w.*?\w)\s*$', 'tokens', 'once', 'lineanchors')); % Find var number numVariables = … str2double(char(regexp(rawfile(simNoVariablesIni(i):simNoPointsIni(i)), '^No\. … Variables:\s*(\d+)\s*$', 'tokens', 'once', 'lineanchors'))); % Find var names simVariablesIni(i) = … simNoPointsIni(i)+11+strfind(rawfile(simNoPointsIni(i)+12:simNoPointsIni(i)+35),' … Variables:'); vartext = rawfile(simVariablesIni(i)+12:simValuesIni(i)-3); variables = regexp(vartext, … '^\s*(?<num>\d+)\s+(?<name>\S+)\s+(?<type>.*?)\s*$', 'names', 'lineanchors'); if numVariables < 1 error('The RAW file:\n\t%s\ncontains no variables.', rawfilename); end if numVariables ~= length(variables) error('Was expecting %d variables, but found %d in RAW file:\n\t%s', … numVariables, length(variables)); end data(i).numVariables = numVariables; % Tanslate vars to valid names in Matlab varnames = {variables.name}; varnames = regexprep(varnames, '(', ''); varnames = regexprep(varnames, ')', ''); varnames = regexprep(varnames, '-', ''); varnames = regexprep(varnames, '#', 'hash'); varnames = regexprep(varnames, '\.', 'dot'); varnames = regexprep(varnames, '\+', 'plus'); varnames = regexprep(varnames, ':', ''); % Remove colons. data(i).varnames = varnames; % Find num points if (strcmp(simType{i},'tran'))||(strcmp(simType{i},'ac')) numPoints = … str2double(char(regexp(rawfile(simNoPointsIni(i):simVariablesIni(i)), '^No\. … Points:\s*(\d+)\s*$', 'tokens', 'once', 'lineanchors'))); if numPoints < 1 145 error('The RAW file:\n\t%s\ncontains no points.', rawfilename); end data(i).numPoints = numPoints; end % next sim if i==nSim nextSim = length(rawfile); else nextSim = simTitleIni(i+1)-3; end % Show vars to save?? if isempty(str2num(outVars{1})) flagShow = 1; ind_oV = 1; else flagShow = str2num(outVars{1}); ind_oV = 2; end switch simType{i} case 'op' % Read data opData = sscanf(rawfile(simValuesIni(i)+10:nextSim),'%f'); data(i).results = opData; if flagShow % Show data text2show = ['Results Operating Point Analysis: ' newline]; numOutVars = length(outVars); for j=ind_oV:numOutVars index = find(strcmp(data(i).varnames,outVars{j})); if isempty(index) warning(['Variable ' outVars{j} ' does not exist.']); else if strcmp(variables(index).type,'current') unit = 'A'; elseif strcmp(variables(index).type,'voltage') unit = 'V'; 146 end text2show =[text2show varnames{index} ' = ' … num2str(opData(index)) ' ' unit newline]; end end fprintf(text2show); end case 'tran' % Read data data_matrix = zeros(numVariables, numPoints); tranData = sscanf(rawfile(simValuesIni(i)+9:nextSim),'%f'); for j=1:numPoints for k=1:numVariables data_matrix(k,j) = tranData((j-1)*(numVariables+1)+1+k); end end data(i).time = data_matrix(1,:); data(i).results = data_matrix(2:end,:); data(i).varnames = data(i).varnames(2:end); % Show data numOutVars = length(outVars); for j=ind_oV:numOutVars index = find(strcmp(data(i).varnames,outVars{j})); if isempty(index) warning(['Variable ' outVars{j} ' does not exist.']); else if strcmp(variables(index+1).type,'current') unit = 'A'; elseif strcmp(variables(index+1).type,'voltage') unit = 'V'; end data(i).outVars(j-ind_oV+1).name = outVars{j}; data(i).outVars(j-ind_oV+1).value = data(i).results(index,:); if flagShow figuremax; plot(data(i).time,data(i).results(index,:)); ylabel([outVars{j} ' [' unit ']']) 147 xlabel('Time [s]'); end end end case 'ac' % Read data rawfileValues = regexprep(rawfile(simValuesIni(i)+9:nextSim),',',' '); mag = zeros(numVariables, numPoints); phase = zeros(numVariables, numPoints); tranData = sscanf(rawfileValues,'%f'); for j=1:numPoints for k=1:numVariables n1 = tranData((j-1)*(2*numVariables+1)+1+2*k-1); n2 = tranData((j-1)*(2*numVariables+1)+1+2*k); cNum = n1+1j*n2; phase(k,j) = angle(cNum)*180/pi-360;%*(n1<1e-3)*(n2>0); mag(k,j) = abs(cNum); end end data(i).f = mag(1,:); data(i).mag = mag(2:end,:); data(i).phase = phase(2:end,:); data(i).varnames = data(i).varnames(2:end); % Show data numOutVars = length(outVars); for j=ind_oV:numOutVars index = find(strcmp(data(i).varnames,outVars{j})); if isempty(index) warning(['Variable ' outVars{j} ' does not exist.']); else if strcmp(variables(index+1).type,'current') unit = 'A'; elseif strcmp(variables(index+1).type,'voltage') unit = 'V'; end data(i).outVars(j-ind_oV+1).name = outVars{j}; data(i).outVars(j-ind_oV+1).mag = data(i).mag(index,:); 148 data(i).outVars(j-ind_oV+1).phase = data(i).phase(index,:); if flagShow figuremax; subplot(2,1,1);semilogx(data(i).f,data(i).mag(index,:)); title(['Bode ' outVars{j}]); ylabel(['Magnitude [' unit '/' unit ']']); subplot(2,1,2);semilogx(data(i).f,data(i).phase(index,:)); ylabel(['Phase [º]']) xlabel('Frequency [Hz]'); end end end otherwise error('Simulation type not supported.'); end end % Remove rawfile delete(rawfilename); % Delete temp netlist % Output Data out = data; end Obtain amplitude and Frequency without Fast Fourier Transform function out = getFA(y,xq,L) flag = 1; ind1 = 1; f = []; A = []; while flag==1 flag1 = 0; flag2 = 0; 149 flag3 = 0; for i=ind1:L if (flag1==0)&&(y(i)>0)&&(flag2==0) flag1 = 1; end if flag1&&(y(i)<0)&&(flag2==0) ini = i; flag2 = 1; end if flag2&&(y(i)>0) flag3 = 1; end if (flag3&&(y(i)<0))||(i==L) fin = i-1; break end end if ((L-(fin+1))<(fin-ini)) flag = 0; else ind1 = fin; end indMax = round(mean(find(max(y(ini:fin))==y(ini:fin)))); indMin = round(mean(find(min(y(ini:fin))==y(ini:fin)))); f = mean([f 0.5/abs(xq(indMax+ini)-xq(indMin+ini))]); ysec = y(ini:fin); A = mean([A max(ysec)]); end out.A = A; out.f = f; end