A systems theory approach to describe dynamic coupling at the cell–electrode interface

In vitro electrophysiology commonly relies on recording techniques, which monitor the electrogenic cell activity through electrodes either placed intracellularly or in close contact with the plasma membrane extracellularly. If extracellular electrophysiology allows for the simultaneous recording and stimulation of multiple sites within a neuronal network,^1,2 the amplitude of the recorded signals, i.e., voltages, is smaller compared to intracellular recordings and sub-threshold potentials are unlikely to be detected.^3,4

The state-of-the-art platforms for extracellular recordings are planar microelectrode arrays (MEAs). Here, the interface between the individual cell and the electrode underneath plays a major role in achieving high quality action potential recordings.^5,6

In the last decade, planar MEAs have been engineered with pseudo-3D micro- and nanoprotrusions recapitulating designs and geometries of living neuronal cells, i.e., dendritic spines, and membrane curvature, which trigger local phagocytosis-like events at the plasma membrane.⁷ These pseudo-3D micro- and nanostructures have been found to decrease the average distance (cleft)^8,9 between the plasma membrane and the electrode surface and increase the coverage area at the interface. Both the optimization of the cleft and the coverage area led to increased amplitude extracellular recordings, ultimately achieving the detection of the cell sub-threshold electrical activity.¹ 

This cell–electrode interface is typically modeled through an equivalent electrical circuit and numerical simulations in which experimental data (i.e., extracellular recordings) are fitted by assigning certain numerical values to the different parameters of the model.¹⁰ 

Among these parameters, the cleft and area coverage are generally considered as average values, thus neglecting any time-dependent variation due to the typical biomechanical reshaping of the plasma membrane and formation of adhesion sites in response to the electrode geometry and material.

Furthermore, these circuital models are strictly related to certain experimental conditions and have limited predictive capability.

In this scenario, a systems theory approach is suitable to describe the continuous changes at the cell–electrode interface by modeling its fundamental parameters with time-dependent functions and their boundary conditions.

Here, we present a cell–electrode interface model that exploits systems theory to model and simulate the dynamical coupling of cells with different types of electrodes (either planar or pseudo-3D). In particular, we investigate the variation of cleft and electrode coverage and their effect on extracellular recordings over time. Finally, we define optimal features and parameters to design pseudo-3D electrodes as active elements of microelectrode arrays.

The cell–electrode interface is modeled considering the physical interaction between neuronal cells and electrodes with planar and pseudo-3D configurations.

Figure 1(a) shows the three main phenomena that occur during an extracellular recording: (1) over the time of an action potential, opening and closing of ion channels result in the flow of ions across the cell membrane (black square), (2) ion barriers between the cell and the electrolyte solution are formed (red square), and (3) cell reshaping and mechanical displacements continuously change the intensity of the coupling over time (green squares). Here, these biological phenomena are modeled by an electrical equivalent circuit [Fig. 1(b)]. In particular, two functions are introduced to characterize the dynamical nature of the coupling: cleft and electrode coverage describing the cell–electrode distance and the portion of the surfaces actually involved in the coupling, respectively.

In this context, the proposed time-varying simulation provides the boundary conditions for recording extracellular potentials over time. Finally, the optimization of the model (dimensioning) suggests the ideal requirements for the electrode design [Fig. 1(c)].

The neuron is modeled by a compartmental approach¹¹ with a partition of the cell into non-junctional (free membrane in the medium and not in contact with the electrode) and junctional (membrane in contact with the electrode) domains. Here, the representation of the plasma membrane, ion channels, and the dynamics of the action potential generation and propagation follow the Hodgkin and Huxley model.¹² This model is described by four differential equations that indicate the voltage variation over time between the inside and the outside of the cell and the activation of sodium and potassium channels.

Here, we consider a simplified approach by defining two differential equations that describe the behavior of the membrane voltage along with the activation of the potassium channels (supplementary material Sec. S1). Given that the sodium and potassium channel dynamics are linearly related,¹³ the differential equations that describe the sodium channel behavior can be neglected,

where

The plasma membrane is represented by a capacitance C_mem, defining the accumulation of charges along the membrane, and a membrane potential V_mem, given by the different accumulation of charges between the intracellular and extracellular domains [Fig. 1(b), black square].

The propagation of an action potential relies on the opening and closing of ion channels, with the resulting flow of ions generating currents across the membrane. Here, the inward currents generated represent the input of the system.

Furthermore, I_k describes the voltage-gated persistent potassium current, which depends on the conductance of the channels (g_k) and the Nernst potential (E_k). I_Na is the sodium current, while I_l describes the leak current. Finally, the output of this system is the membrane voltage (supplementary material Sec. S1, Fig. S1).

The second part of model is the equivalent circuit of the junctional membrane, which is the portion of the cell in contact with the electrode [Fig. 1(b), black square]. Here, the plasma membrane is represented by a capacitor (C_j) and a resistor (R_j). Assuming that (1) voltage-gated channels have not been included¹¹ since their presence would not influence the propagation of the action potential and (2) no current is generated in the proximity of the electrode, the initiation of an action potential generates equipotential conditions inside and outside the cell.

Finally, the non-junctional and junctional components are electrically connected through the axoplasmic resistance that is representative of the ion flow inside the cell volume propagating the voltage spikes across the membrane.¹⁴ 

At the interface between the cell and the electrode, an interspace (cleft) is formed and filled with an extracellular medium.¹⁵ The presence of this electrolyte solution has three effects: (1) accumulation of charges at the cell–electrolyte solution interface, (2) electrical power dissipation within the cleft, and (3) accumulation of charges at the electrolyte solution–electrode interface. The accumulation of charges both at the cell–electrolyte interface and at the electrolyte solution–electrode interface is characterized by a double barrier of ions, called electrical double layer (EDL), where charges are accumulating both at the cell and electrode surfaces.¹⁶ Here, the cell–electrolyte solution and the electrolyte solution–electrode interface are represented by a capacitor and a resistor connected in parallel [Fig. 1(b), red square, R_h–C_h and R_e–C_e]. The capacitance models the accumulation of charges, due to the presence of the ions along the surfaces of the cell and the electrode, while the resistance describes the opposition to the ion flow across the EDL.

Both the average cleft width and the extracellular medium presence directly affect the amplitude of the recorded signals, influencing the local field generated by the cell electrical activity. Here, the power loss of the recorded signal is modeled as a (sealing) resistance R_seal¹⁷ [Fig. 1(b), green square].

The measurement stage is modeled as a capacitor and a resistor connected in parallel (R_m and C_m, respectively), while the electrode resistance contribution can be neglected as previously shown.^1,18

The proposed electrical equivalent circuit is mainly dominated by a resistive–capacitive behavior. Thus, the voltages of the capacitors (i.e., junctional membrane, EDL layers, and measured voltage) are dynamic variables. The mathematical model describes the behavior of these voltages as a function of time by using differential equations.

Applying Kirchhoff’s circuit laws, the model is defined as follows:

Here, each equation describes the total current flowing through the capacitors. To study this model in a state-space form, we define the state vector, i.e., the set of variables needed to describe the dynamics of the system, as follows:

V_mem is the membrane voltage described in (1), modeling the propagation of an action potential across the cell. This represents the input of the “interface system” and is referred to as u(t)—following the standard notation of systems theory. The output of the system, defined as the set of physical quantities to be observed, is referred to as y(t).

The state-space model exploits matrices to express the set of equations that describes the system under study. The time-derivative of the state vector is $ẋ(t)$⁠, and here, four matrices are introduced: (1) matrix A relates the state vector to its derivative (i.e., its evolution over time), (2) matrix B describes how the input is directly affecting the evolution of the state vector, (3) matrix C links the state vector to the output of the system, and (4) matrix D describes the relationship between the input and the output. The state-space model is eventually described as follows:

Notably, matrix A(t) is the dynamical matrix of the system and includes the variation of the sealing resistance, allowing the simulation of a dynamical coupling.

Further details on both the mathematical model and state-space representation are discussed in supplementary material Sec. S2.

Here, numerical values are assigned to the characteristic parameters of the model (i.e., geometrical dimensions of the electrodes, EDL equivalent parameters, etc.). For instance, when planar electrodes are adopted at the recording site, their radii and resulting surface are considered (Table I).

In the case of pseudo-3D electrodes, mushroom- and pillar-like shapes have been studied. In mushroom-shaped microelectrodes, the effective electrode surface is the sum of a cylinder later surface (stalk) with radius in the range of 0.3 µm–2 µm and height in the range of 0.5 µm–2 µm and a semi-sphere surface (cap) with radius in the range of 0.6 µm–2.5 μm.¹⁹ In the case of pillar-like electrodes, a cylinder with radius in the range of 0.15 µm–1.5 µm and height in the range of 0.3 µm–3 μm has been considered.

The resistance of the EDL at the cell–electrolyte solution interface (R_h) is obtained from the experimental measurements previously shown.²⁰ 

The resistance and the capacitance of the measurement stage are chosen as follows: R_meas = 10 000 GΩ and C_meas = 8 pF.¹⁸ These values are typical amplifier input values, accounting for an impedance at the input stage of 20 MΩ at 1 kHz.

The neuronal cell is modeled by a sphere (soma) of 30 µm radius and a cylinder (axon) of 3 µm radius and 150 µm length. This surface is divided into non-junctional and junctional areas, where the latter strongly depends on the electrode surface in contact with the cell.

Importantly, in the case of planar electrodes, given that cells change their physical configuration and adhesion processes during the coupling, the junctional surface might be larger than the electrode area. For this reason, the junctional surface is chosen as the double of the electrode surface.

Similarly, in the case of pseudo-3D electrodes, cells gain intimate adhesion to the electrode surface, and therefore, the resulting junctional surface might be evaluated as 30% of the stalk surface plus cap and rim surfaces for mushroom-shaped microelectrodes, and 30% of the stalk surface plus the top side in the pillar case as previously shown.¹⁹ 

Moreover, the numerical value of the junctional surface is used in the computation of both EDL parameters—R_h and C_h—and junctional membrane parameters C_j and R_j (Table I), while the non-junctional surface area is used to compute Hodgkin and Huxley model parameters as shown in supplementary material Sec. S1.

Furthermore, the modulation of the junctional surface allows for the simulation of different coupling conditions, such as the coupling occurring at the axonal domain, as shown in supplementary material Sec. S1, Fig. S2. Finally, the flow of ions during the propagation of an action potential can be modeled as a current flowing through a (axoplasmic) resistance, which is proportional to the total volume of the non-junctional cell.²¹ 

Tables I and II summarize the nominal values assigned to the aforementioned parameters in the case of individual planar and 3D electrodes, respectively. More details on the dimensioning are reported in supplementary material Sec. S3.

As mentioned above, the sealing resistance is a fundamental parameter in the cell–electrode interface, regulating the current leakage between the neuron and the electrode. It is influenced by both the cleft width and the presence of the electrolyte solution and its value directly affects the amplitude of the extracellular recordings.

For planar electrodes, the sealing resistance is estimated as follows:²¹ 

where ρ_s is the resistivity of the electrolyte solution, d is the average cleft distance, and δ is a coverage percentage, which models the fraction of electrode surface actually in contact with the neuron.

In pseudo-3D electrodes (both pillar and mushroom-shaped), the planar contribution to the computation of the sealing resistance is neglected, assuming that they are electrically isolated. In the pillar-like electrode, the sealing resistance is the sum of two contributions given by the adhesion areas on the pillar stalk and top area,²² 

where δ_stalk and δ_top are coverage percentages for the stalk and the top of the pillar-shaped electrode, respectively. The stalk coverage coefficient may change, while the coverage of the top surface is fixed to 0.3.¹⁹ 

In the mushroom-shaped electrode, the sealing resistance takes into account three different contributions,¹⁹ 

Coverage percentages δ_stalk and δ_rim are introduced to describe the percentage of surface involved in the computation of the sealing resistance.

Table III summarizes the nominal sealing resistances estimated in this model, considering a group of four vertical structures in the case of pseudo-3D electrodes.

Numerical simulation, computation, and curve fittings were performed with Matlab and Simulink. Details on simulations and general numerical procedures are shown in supplementary material Secs. S5 and S8.

A significant membrane reshaping may occur when cells are in contact with the electrode underneath, and this might affect extracellular recordings of action potentials over time. This reshaping causes a point-by-point change of the cleft forming between the plasma membrane and the electrode surface. Therefore, we consider a cleft function d(t), which can be arbitrary defined in order to mimic a vertical displacement of the cell approaching or moving away from the electrode. The cleft function is then used to compute and predict the sealing resistance of the model.

Figure 2(a) shows how the sealing resistance changes as a function of the cleft distance when a cell adheres on planar, pillar, and mushroom-shaped electrodes (blue, red, and yellow lines, respectively). The cleft distance and sealing resistance are inversely correlated [Eqs. (5)–(7)], resulting in high resistance values when the cleft distance tends to zero, while showing a (almost) steady value for larger cleft distances (more than 80 nm).

Figure 2(b) shows a time-varying simulation of an extracellular recording in which the cleft distance linearly increases over time. Simulated extracellular recordings for the planar electrode are less sensitive to cleft variations compared to pillar and mushroom-shaped electrodes, resulting in a much slower decay of the amplitude of the simulated signal. The recordings shown in Fig. 2(b) are computed by assuming an increase of the cleft distance; however, the model is yet valid for any chosen behavior of the cleft distance and surface coverage functions (supplementary material Sec. S4, Figs. S5–S7).

Supplementary material Sec. S5 also provides details on the computation and numerical simulation procedure.

Here, the simulations are computed on a representative time frame of 120 ms; however, the systems theory approach allows for simulations extended to any relevant time frame. This is fundamental to span from short term cell–electrode coupling phenomena to achieve electrophysiological recordings, which are typically carried over months. More details can be found in supplementary material, Fig. S10.

Cells coupled to electrodes may also move from side to side with respect to the recording electrode or may bend. These micro-movements might reduce the effective surface where the cell and electrode are physically coupled. For this reason, surface coverage percentages were defined in Eqs. (5)–(7). As for the cleft distance, a static parameter is not representative of the continuous reshaping of the plasma membrane; therefore, here, the function coverage(t) is introduced to ideally model the variation over time of the physical coupling between the cell and the electrode.

Figure 3(a) depicts the sealing resistance as a function of the coverage percentage, considering planar, pillar, and mushroom-shaped electrodes. Here, the relationship between the sealing resistance and coverage is linear [see also Eqs. (5)–(7)]. By comparing Figs. 2(a) and 3(a), it is possible to note that the cleft average distance influences more the quality of the recordings compared to the coverage percentage. The cleft distance, in fact, induces a hyperbolic behavior of the sealing resistance that is significantly reduced when the cleft average distance increases.

The coverage function, instead, is related to a linear variation, thus inducing a proportional dampening action of the electrical signals at low coverage values. This consideration is further validated by simulated recordings [Fig. 3(b)] as the recorded signal amplitude linearly increases with the coverage for all electrodes. More details can be found in supplementary material Sec. S4, Fig. S7.

Previous analysis showed how the variation of cleft or coverage functions individually could influence the resulting extracellular recording. However, the dynamic coupling process is effectively described by the combination of both vertical and lateral membrane reshaping. Here, we evaluated the simultaneous effect of the cleft and coverage functions to determine the sealing resistance variation on extracellular recordings over time [Fig. 4(a)].

The sealing resistance, in the case of planar electrodes, is strongly dependent on the coverage: high values of coverage induce the highest values of sealing resistance.

In contrast, in the pseudo-3D electrode case, a small cleft and low coverage [Fig. 4(a), blue lines] may lead to higher values of the sealing resistance. Supplementary material Sec. S7 (Figs. S11–S13) shows the sealing resistance representation as a function of both the cleft and the coverage.

Time-varying simulation of extracellular recordings, in which both cleft and coverage functions change over time, shows the strict correlation of these two parameters with the amplitude of the recorded signals. More specifically, increasing the coverage and/or decreasing the cleft distance will generate an increase in the amplitude of simulated recordings, whereas decreasing the coverage and/or increasing the cleft distance will decrease the amplitude of such recordings, probably due to cell–electrode decoupling [Fig. 4(b) and supplementary material Sec. S7, Figs. S11–S13].

Finally, time-varying simulation is exploited to predict a detection limit of each electrode type.

Considering voltage thresholds of 400 μV, 200 μV, and 50 μV (for planar, pillar, and mushroom-shaped electrodes, respectively), we evaluated the conditions of the cleft distance and coverage functions for the different electrode types necessary to achieve the desired extracellular recording [Fig. 4(c)].

In addition, Fig. 4(c) depicts how planar electrodes rely on higher values of coupling parameters, compared to pseudo-3D electrodes. For instance, considering a coverage of 40%, only pseudo-3D electrodes can record extracellular potentials of the desired voltage, whereas planar electrodes would need a coverage of at least 50%. If the model is fitted on real recordings, the computation of the detection limits may provide information and predictions on actual experiments.

It is important to highlight that these results prove the validity of the mathematical modeling that is used to predict the electrode behavior. An a priori fitting of the model to actual recordings allows for a realistic computation of the predicted behaviors.

Supplementary material Sec. S9, Fig. S18 illustrates the simulated electrical potential when the detection limits of the electrodes are not satisfied (i.e., in the case of high cleft distances).

Pseudo-3D electrodes actively change the mechanical membrane response, promoting an engulfment-like mechanism at the cell–electrode interface that affects the recorded extracellular signal. This mechanism is correlated with the membrane bending as response to the pseudo-3D electrode shape. The main geometrical parameter that has an impact on the membrane wrapping is the aspect ratio of the structures.

The engulfment mechanism is a variable phenomenon depending on the cell type; however, it is mainly described by the combination of the two coupling functions.

Therefore, a precise design of the 3D geometry is indeed essential in order to optimize extracellular recordings. Figure 5(a) schematically shows possible membrane responses to different microprotrusions at different aspect ratios.^22,23 Here, the engulfment percentage is defined as the portion of the electrode surface wrapped by the plasma membrane. Figure 5(b) shows the percentage of the engulfment of a pillar electrode and the engulfment of the stalk and the rim in a mushroom-shaped electrode, derived from previous experimental work.²² 

Moreover, the total junctional surfaces and sealing resistances are computed as a function of the geometry prior to the numerical optimization process (Table IV and supplementary material Sec. S8, Figs. S14–S17).

For the pillar electrode, the optimization of the aspect ratio increases the contact surface. The increased contact surface enhances the sealing resistance of 7% and consequently improves the overall recorded signal of about 70%. Likewise, the optimization of the aspect ratio for mushroom-shaped electrodes will cause a 10% increase of the sealing resistance and a 37% increase of the overall recorded signal.

The numerical values obtained from this optimization process can be exploited as guidelines for the design and fabrication process of a pseudo-3D electrode.

The cell–electrode interface is typically modeled by an electrical equivalent circuit. However, such circuital models do not take into account the highly dynamic mechanical behavior of cells and can partially simulate the electrical properties of neuronal cells coupled to microelectrodes.

In this work, we introduced a mathematical model and a systems theory approach to numerical simulations, in order to characterize the cell–electrode interface through the vertical (cleft function) and lateral (electrode coverage function) displacements of the membrane over time. Unlike conventional circuital approaches, the mathematical model involves both electrical and mechanical conditions of cells over time.

We demonstrated the correlation between the dynamic coupling and extracellular recordings to ultimately define the detection limits of planar and pseudo-3D electrodes, based on experimental data. The implementation of a bio-mechanical model of the membrane could represent an improvement for this approach.

Finally, the large variability of extracellular recordings, induced by the variation over time of the coupling, is analyzed.

The final optimization of the electrode geometrical parameters and the prediction of the potential coupling might be exploited for future microelectrode designs prior to experimental validation.

This time domain model lays the ground work to further investigate the more refined cleft-coverage dynamics of cells in light of the emerging nano-patterned electrodes.

See the supplementary material for details on numerical fittings and mathematical procedures. Additional simulations are shown in order to clarify the findings of this article.

We thank Csaba Forro, Valeria Criscuolo, and Claudia Lubrano for comments and advice on the manuscript.

The data that support the findings of this study are available within the article and its supplementary material.