State-space modeling for dynamic response of graphene FET biosensors

Since the introduction of ion-sensitive field effect transistors (ISFETs) in 1970,¹⁾ FET-based biosensors have attracted significant interest owing to their promising potential for a wide range of applications.^2,3) In particular, graphene films offer an ideal sensing platform owing to their high electron/hole mobilities⁴⁾ and 2D nature,⁵⁾ and thus graphene field effect transistor (G-FET) biosensors have been used to detect ions,⁶⁾ biomolecules,⁷⁾ and bacteria.⁸⁾ However, a baseline drift is observed in the response of FET-based sensors,^9,10) particularly in the case of G-FET biosensors because of their high sensitivity. The baseline drift makes it difficult to accurately estimate the concentration of target molecules. Although various approaches have been proposed to overcome this issue,^11,12) the mechanism of baseline drift remains unclear. The uncertainties in the measurement system, such as the temperature, ion concentration, and pH, may also affect the sensor signal. Therefore, it remains challenging to compensate for the signal using hardware and/or describe the signal using simple models.

We have previously proposed simple state-space models (SSMs) to describe the time-series data of a G-FET biosensor.¹³⁾ Herein, we present more sophisticated SSMs and provide detailed discussions. State-space modeling is a framework established to understand stochastic and deterministic dynamical systems, referred to as states, which are observed through a stochastic process.¹⁴⁾ SSMs are widely used in time-series analyzes, e.g. in Apollo and Polaris aerospace programs¹⁵⁾ and for studying animal movements.¹⁴⁾ Least-squares methods¹⁶⁾ are unsuitable for the state and parameter estimation of SSMs, because they systematically under- or overestimate unknown model parameters, as there is a serial correlation between successive observations.¹⁷⁾ The Kalman filter (KF) algorithm is commonly used for the parameter estimation of SSMs.¹⁸⁾ However, in the KF algorithm, the state and parameter are represented by a Gaussian probability distribution function. This makes it unsuitable for nonlinear dynamical systems. Therefore, we used the Markov chain Monte Carlo (MCMC) algorithms for state and parameter estimation.¹⁹⁾ The MCMC algorithms provide a general methodology that can be applied to nonlinear and non-Gaussian state models, because they allow sampling from an arbitrary posterior distribution. After estimating the state and parameters through the MCMC methods, we divided the time-series data of the G-FET biosensor into sensor response to target molecules and baseline drift through the proposed SSMs. In addition, we built competitive SSMs to determine the one that best describes the obtained dynamic response of the G-FET biosensor.

Figure 1(a) shows a schematic of the measurement system including a G-FET sensor. A graphene film was first grown on a Cu foil by chemical vapor deposition and then transferred onto a Si/SiO₂ substrate. A source/drain electrode was formed with 10 nm Ti and 90 nm Au. Finally, a graphene channel was formed by O₂ plasma etching. The obtained G-FET was immersed in 15 mM D-PBS (Nacalai Tesque, Inc.). A bias voltage (V_DS) of 0.1 V was applied, and a top-gate voltage (V_GS) was applied through the solution using a Ag/AgCl electrode. While sweeping V_GS, the drain current (I_DS) was measured using a semiconductor parameter analyzer (Keysight Technologies, B1500A). Figure 1(b) shows a typical I_DS–V_GS plot. From the I_DS–V_GS plot, the charge neutral point (CNP), which is V_GS at minimum I_DS in the plot, was calculated. Figure 1(c) shows the time-series data of the CNP. While monitoring the CNP, bovine serum albumin (BSA) molecules (Sigma-Aldrich), as a sensing target, were intermittently exposed to the sensor, with varying concentrations of 50 pM, 500 pM, 5 nM, 50 nM, 500 nM, 5 μM, and 10 μM, in D-PBS solution [Fig. 1(d)]. The BSA molecules were physisorbed on graphene. In the dynamic response of the sensor, shown in Fig. 1(c), the CNP varies with the BSA concentration. In addition, there is an apparent baseline drift, as the CNP sequentially varies even when there is no BSA in the solution.

Zoom In Zoom Out Reset image size

3.1. Results

To model the obtained time-series data in Fig. 1(c), we propose an SSM, shown in Fig. 2(a), described by Eqs. (M 1.1)–(M 1.6). We define the model as Model 1

$$\begin{eqnarray}&&{{y}}_{t}={{x}}_{t}+{{q}}_{t}+{\varepsilon }_{t},\,\,{\varepsilon }_{t}\sim N\left(0,{\sigma }_{\varepsilon }\right),\end{eqnarray} \tag{ M 1.1 }$$

$$\begin{eqnarray}&&{{x}}_{t}-{{x}}_{t-1}={{x}}_{t-1}-{{x}}_{t-2}+{{w}}_{t},\,\,{{w}}_{t}\sim N\left(0,{\sigma }_{{w}}\right),\end{eqnarray} \tag{ M 1.2 }$$

$$\begin{eqnarray}&&{{\rm{q}}}_{t}=\displaystyle \sum _{i=1}^{N}{{\rm{\Delta }}q}_{i}\end{eqnarray} \tag{ M 1.3 }$$

$$\begin{eqnarray}&&{{\rm{\Delta }}q}_{i}=\left\{\begin{array}{l}0\,\,\,\,\,\,\,\,\,\,\,\,\,\,\left(t\lt {t}_{i}\right)\\ \left({Q}_{i}-{Q}_{i-1}\right)\left\{1-\exp \left(-t-{t}_{i}/{\tau }_{i}\right)\right\}\,\,\,\,\left(t\geqslant {t}_{i}\right)\end{array}\right.,\end{eqnarray} \tag{ M 1.4 }$$

$$\begin{eqnarray}&&{{Q}}_{i}=\displaystyle \frac{a{{c}}_{i}}{{{K}}_{D}+{{c}}_{i}},\end{eqnarray} \tag{ M 1.5 }$$

$$\begin{eqnarray}&&{\tau }_{i}=A\,\mathrm{log}\left({{\rm{c}}}_{i}\right)+B+{\omega }_{i},\,{\omega }_{i}\sim N\left(0,{\sigma }_{\tau }\right).\end{eqnarray} \tag{ M 1.6 }$$

Zoom In Zoom Out Reset image size

According to Eq. (M 1.1), which is called the observation equation, the observed CNP (y_t) is a summation of the baseline (x_t), regression component regarding BSA (q_t), and observation noise [ε_t ∼ N (0, σ_ε)]. The subscript t (= 1, 2, 3, ..., 104) corresponds to the data index. In Eq. (M 1.2), which is called the state equation, the unobserved baseline (x_t) is assumed to follow a quadratic trend with the system noise [w_t ∼ N (0, σ_w)]. The regression component (q_t) is described based on two concepts. The first concept is that the signal reaches a plateau defined by the Langmuir model in Eq. (M 1.5), where K_D, c_i, and a are the dissociation constant, BSA concentration, and a coefficient, respectively. The subscript i (=0, 1, 2, ..., 7) is the sample index, which increases when BSA is introduced to the sensor. The assumption, represented by Eq. (M 1.5), agrees with experimental data,^6,20) confirming that the sensor response to BSA follows the Langmuir equation. The shift direction of the CNP against the adsorption of BSA is debatable. Some studies showed that the CNP shifts in the positive direction,⁶⁾ whereas others showed that it shifts in the negative direction.^21,22) This is because the detection mechanism involves electrostatic gating effect,²³⁾ charge doping,²⁴⁾ and charged-impurity scattering.²⁵⁾ The difference is presumably due to the surface condition of graphene, which is influenced by the manufacturing process. According to our data, shown in Fig. 1(c), the CNP seems to shift in the negative direction. The second concept is that the signal varies following an exponential decay after the BSA is introduced [Eq. (M 1.4)]. According to previous studies, the response time of FET-based biosensors is in the order of seconds to minutes, even though the target concentration ranges from the order of fM to the order of μM.^26,27) This implies that the decay time τ_i is not directly proportional to the concentration (c_i). Therefore, in Eq. (M 1.6), we assume that the mean of τ_i is proportional to the logarithm of c_i, not directly proportional to c_i.

We estimate the states and parameters of Model 1 to fit the dynamic response shown in Fig. 1(c). As the regression component (q_t) in our SSMs is highly nonlinear, the KF algorithms cannot be applied to the subsequent parameter estimation. Therefore, we used the MCMC methods¹⁹⁾ implemented in Python and Stan's probabilistic programming languages²⁸⁾ to estimate the state and parameters of Eqs. (M 1.1)–(M 1.6). In the computation, four chains were run with 10 000 iterations for each chain. The first 5000 iterations were discarded as burn-in, and the last 5000 iterations in each chain were used for subsequent inferences. Figure 2(b) shows the trace plots for K_D. The posterior distribution is well sampled, as there is no serial correlation, and the chains explored the sample space many times. The other model parameters have similar trace plots (some are shown in Fig. S1, available online at stacks.iop.org/JJAP/59/SGGH04/mmedia). The Rhat value, which is an indicator of model convergence,²⁹⁾ was lower than 1.1 for all the parameters (Table I and Table SI). This result indicates that the model reaches convergence. The histogram of the trace plots, shown in Fig. 2(b), corresponds to the posterior distribution of the dissociation constant K_D [Fig. 2(c)]. The posterior mean and the 95% credible interval for K_D were calculated to be 151 nM and (85 nM, 257 nM), respectively. These values are comparable to those reported previously.⁶⁾ Table I lists some representative parameters (see other parameters in Table SI). Using the proposed SSM and estimated parameters, we divided the time-series data, shown in Fig. 1(c), into the baseline [x_t, Fig. 3(a)] and a signal related to BSA [q_t, Fig. 3(b)]. Figure 3(c) shows the sensor response against BSA as a function of the concentration. It should be noted that the effect of baseline drift is excluded in Fig. 3(c), and thus it is indicated that the dynamic response in Fig. 1(c), in which the BSA concentration ranges from 0 to 5 nM, was mainly attributed to the baseline drift.

Table I.  Summary of representative parameters for the SSM estimated using the MCMC method. The posterior mean, the 95% credible interval, and Rhat are shown for K_D, a, σ_ε, and τ₆. τ₆ corresponds to the decay time when the BSA concentration increases from 500 nM to 5 μM.

	Mean	2.5%	97.5%	Rhat
KD [nM]	151	85	257	1.0
a [mV]	−81	−120	−53	1.0
σε [mV]	2.2	1.9	2.6	1.0
τ6 [min]	6.35	1.46	13.83	1.0

Zoom In Zoom Out Reset image size

3.2. Discussion

We developed SSMs to describe the dynamic response of a G-FET biosensor against BSA. The parameters of the models were estimated using the MCMC methods. The SSMs effectively fitted the dynamic response. Using the estimated parameters, we divided the time-series data into BSA response and baseline drift. Although the true baseline or states cannot be observed directly, state-space modeling helps estimate the probabilities of the states and parameters by explicitly including the observation equation. Furthermore, it allows incorporating any established theories and/or empirical insights into the models, making it flexible. Although the true state equation is unknown, state-space modeling can provide an insight into the dynamic data by generating a series of competing models and comparing them using model selection methods based on criteria such as the WBIC. The proposed method can be applied not only to G-FET biosensors but also to FET-based sensors for an accurate analysis of the sensor response. Therefore, this work paves the way for further applications based on FET-based sensors.

This work was supported by JST CREST Grant Number JPMJCR15F4, Japan.