A.
Toral-Lopez
*a,
D. B.
Kokh
b,
E. G.
Marin
a,
R. C.
Wade
bcd and
A.
Godoy
a
aDepartamento de Electrónica y Tecnología de Computadores, Facultad de Ciencias, Universidad de Granada, Spain. E-mail: atoral@ugr.es; agodoy@ugr.es
bMolecular and Cellular Modeling Group, Heidelberg Institute for Theoretical Studies, Schloss-Wolfsbrunnenweg 35, 69118 Heidelberg, Germany
cCenter for Molecular Biology (ZMBH), DKFZ-ZMBH Alliance, Heidelberg University, Im Neuenheimer Feld 282, 69120 Heidelberg, Germany
dInterdisciplinary Center for Scientific Computing (IWR), Heidelberg University, Im Neuenheimer Feld 205, Heidelberg, Germany
First published on 17th June 2022
Biological Field-Effect Transistors (BioFETs) have already demonstrated enormous potential for detecting minute amounts of ions and molecules. The use of two-dimensional (2D) materials has been shown to boost their performance and to enable the design of new applications. This combination deserves special interest in the current pandemic caused by the SARS-CoV-2 virus which demands fast, reliable and cheap detection methods. However, in spite of the experimental advances, there is a lack of a comprehensive and in-depth computational approach to capture the mechanisms underlying the sensor behaviour. Here, we present a multiscale platform that combines detailed atomic models of the molecules with mesoscopic device-level simulations. The fine-level description exploited in this approach accounts for the charge distribution of the receptor, its reconfiguration when the target binds to it, and the consequences in terms of sensitivity on the transduction mechanism. The results encourage the further exploration of improved sensor designs and 2D materials combined with diverse receptors selected to achieve the desired specificity.
BioFETs take advantage of the same operation principle of well-established Metal-Oxide-Semiconductor FETs (MOSFETs) to detect ions or biomolecules immersed in a solution. In a regular MOSFET, the charge in the semiconductor channel is modulated by the bias applied to the metal gate over it, which is electrostatically coupled through an oxide. Thus, the gate tunes the channel conductivity, electrically joining/isolating the two contacts located at its edges, i.e. source and drain, and consequently controlling the current flow between both ends. In a BioFET, instead, the gate metal is substituted by a liquid solution containing the ions or molecules to be detected, while the insulator in direct contact with the electrolyte is functionalized to enhance the sensitivity and selectivity of the sensor. Then, the charged analytes attached onto the surface induce a change in the semiconductor surface potential that is transduced by a variation of the source-to-drain current. Ultimately, any alteration in the charge adhered to this layer could be detected as a modification in the device output current. Although this operation principle dates back to the 70s,8 huge research efforts in the last decade have significantly improved its performance, and this along with a continuous price reduction has facilitated BioFETs' widespread use.
These efforts have recently been boosted by the surge in interest in graphene and two-dimensional materials (2DMs), which has had a profound impact on sensor design, showing enormous potential to boost sensor performance.9 However, in order to leverage the combination of both technologies and in particular their use for a fast, reliable and cheap detection of ions and molecules, numerous challenges still need to be overcome. While from the fabrication and experimental perspective, a wide variety of results are readily available,10–14 there is a noticeable lack of theoretical understanding concerning the mechanisms and principles ruling their performance. There are several approaches based on both ad hoc and TCAD solutions, but these are focused on nanowire BioFETs (NW-BioFETs)15–17 or planar devices relying on bulk semiconductors.18–20 Only a few are meant to explore 2DM-based BioFETs.21–23 In addition, most of them are based on coarse approximations, such as: (i) treating the sensing layer as a uniformly charged layer or discrete box-shaped charged blocks for the biomolecules, or (ii) analytical or simplified models are employed for the electrolyte. The main reason behind this theoretical gap is the lack of a multidisciplinary and multiscale description of the device operation, necessary to achieve a holistic picture able to consistently capture features ranging from molecular interactions at the solid–liquid interface to the extracted characteristics. This is of particular importance to comprehend the subtleties associated to the sensing target, and predict in an accurate way the sensitivity of the BioFET for the specific objective.
In this context, this work presents a novel multiscale scheme that combines a precise modelling of the receptor and target molecules at the sensor surface and a sophisticated description of the semiconductor device and the liquid electrolyte. Our approach pursues the integration of information from detailed computations of protein electrostatic properties at the atomic level into the device level numerical simulations. This approach provides a potentially powerful tool to understand in detail the operation principle and unveil the optimization keys of Graphene based BioFETs (GBioFETs) for the early detection of SARS-CoV-2.
The rest of the paper is organized as follows. First, we introduce the multiscale methodology exploited to simulate the device with special attention to the sensing layer. Next, we focus on the atomistic modelling of the molecules involved in the detection of SARS-CoV-2. Finally, the results obtained for the GBioFET are presented and analysed and the main conclusions are drawn.
In particular, in the channel, the continuity equation is solved in the diffusive regime assuming a common pseudo-Fermi level for electrons and holes due to the short recombination times measured in graphene.24–26 The drift-diffusion scheme provides the most appropriate description according to the dimensions of state-of-the-art experimental 2D BioFETs, where carrier transport properties are dominated by scattering processes. The source-to-drain current IDS can be expressed as:
IDS = q(nμn + pμp)∇VEF | (1) |
At the electrolyte, we solve the modified Poisson–Boltzmann equation that takes into account the finite size of the ions that are present in the electrolyte region, limiting the maximum ion concentration.28 Using the standard Poisson–Boltzmann equation, on the contrary, can give rise to non-physical values, particularly for large ions. We complete this approach with a sophisticated model for the interactions between the solution and the device surface that encompasses a position-dependent water permittivity εw near the solid–liquid interface and Potentials of Mean Force (PMFs)29,30 as well as the electrolyte complex reactions (see ESI†). The position-dependent εw profile is included in the solution of the Poisson–Boltzmann equation, while the PMF profiles are integrated in the calculation of the concentration of the i-th ion type:
(2) |
In addition to the ions contained in the solution, the electrolyte–device interface is charged with the receptor molecules attached to the device surface. Each receptor can be in one of two states: idle, i.e. no target molecule is attached to them, or activated, i.e. the target molecule is captured. In both cases the molecule charge is treated taking into consideration its shape and spatial distribution, employing an atomic level description of the molecular structure. This charge is later accurately mapped at the biosensor level. The model of the complex electrolyte combined with small charged receptor molecules has been evaluated against experimental data on the surface potential for a histidine–aspartic acid peptide and Green Fluorescent Protein taken from ref. 32 showing an excellent capability to reproduce them as discussed in the ESI 3.† The integration of the molecules shape, as here described, improves the possibilities to reproduce the experimental selectivity, reducing the equivalence between similarly charged molecules at the device computational level.
In the specific case addressed here, i.e. the detection of SARS-CoV-2, prior to the extraction of the molecule shape and charge distribution in the idle (ρidle) and activated (ρact) states, it is worth recalling the infection mechanism of this virus so to define the receptor-target pair of the biosensor. The structure of this virus is determined by a lipid membrane encasing its genetic load. The most remarkable element in this membrane is a set of spike proteins that surround the whole structure, and give rise to the characteristic appearance of the virus, a circular body enclosed by a halo that looks like a corona. These spike proteins, as indicated in Fig. 1, are composed of two subunits, named S1 and S2, respectively. The former is meant to bind the virus to the cell surface while the latter enables the virus core envelope to join the cell membrane releasing its genetic load.33
Fig. 1 (a) Three-dimensional structure of the Spike glycoprotein head in its closed state (PDBe ID:6vxx34) and open state (PDBe ID:6vyb34). It is composed by two subunits named S1 and S2. In the open state the S1 Receptor Binding Domain (S1RBD) is exposed so that the virus can be attached to the cell surface. (b) Structure of the human Angiotensin Converting Enzyme 2 (ACE2) (PDBe ID:1o86 (ref. 35)) used as attachment point by SARS-CoV-2; and its complex with the S1RBD (PDBe ID:6m0j36). |
The main entry point to the cell used by SARS-CoV-2 is the Angiotensin Converting Enzyme 2 (ACE2),37–39 (see Fig. 1). The ACE2 molecule interacts with the S1 Receptor Binding Domain (S1RBD) so that the virus attaches to the cell surface. After this capturing process, the S1 unit is cleaved from the main body of the spike protein leaving the S2 unit exposed, making possible the fusion of the cell and virus membranes.33 According to the aforementioned infection mechanism, here we focus on the ACE2-S1RBD pair, where the ACE2 acts as receptor molecule.
To achieve an accurate description of the electrostatic properties of the ACE2 and ACE2-S1RBD complex, they are modelled in atomic detail. To do so, we employed structures with PDBe ID: 1o8635 and PDBe ID: 6m0j36 for the ACE2 and ACE2-S1RBD complex respectively. Glycans (NAG) were removed from the structures and missing loops were built using MODBASE web site.40 One Zn2+ and one or two (complex and ACE2, respectively) Cl− ions were retained in calculations. The proteins were protonated at pH 7.4 and partial charges at each atom were generated using PDB2PQR server41 (ions were omitted at this step). The total charge of the ACE2 protein and the complex is −12e and −22e, respectively. Zn2+ and Cl− were added back to the structures.
The resulting charge distribution of the unbound (idle state) ACE2 receptor (ρidle) and the bound (activated state) ACE2-S1RBD pair (ρact) were then mapped in a spatial grid of 0.5 nm × 0.5 nm (Fig. 2) and included in the 2D biosensor simulation (see Fig. S3 in ESI†). These charge profiles were replicated along the device in accordance with the state and positions defined for the receptor molecules.
Fig. 2 The charge distribution of apo-ACE2 (a) and the ACE2-S1RBD (b) complex viewed from the projection plane. Blue dots indicate positively charged elements while red dots indicate negatively charged elements in the molecule structure. Figures (c) and (d) shows the final result of the projection of the 3D charge distributions (see Fig. S.3 in ESI†) prior to their integration in the simulation of the GBioFET device. |
Thus, we test the sensitivity of the GBioFET to the presence of ACE2-S1RBD pairs by determining through self-consistent simulations the transfer response of the device (IDS − VFG) assuming different percentages of activated receptors, α (Fig. 4a).
The device response can be qualitatively split into three regions: (i) the gate bias around the point of minimum conductivity, i.e. the Dirac voltage VDirac, (ii) the p-type branch, corresponding to negative gate biases, i.e. to the left of VDirac, and (iii) the n-type branch, corresponding to positive gate biases, i.e. to the right of VDirac. As can be observed at VDirac, the channel has the lowest conductivity and the GBioFET does not show noticeable changes when the S1RBD molecules are attached to the ACE2 receptors, indicating that this bias region is poorly responsive to detect the presence of the virus. Concerning the p-type branch, a slight modification in the output current is observed as α varies: the transfer response spreads out as the gate bias gets more negative, but within a limited range. The n-type branch reveals itself as the most sensitive operation region: the magnitude of the output current is lower than the p-type branch, but it exhibits a higher sensitivity to changes in α.
More importantly, the transfer characteristic manifests a counter-intuitive behaviour with the activation percentage. Both the ACE2 receptor and the S1RBD are negatively charged. Thus, when the receptors are activated (i.e. the ACE2-S1RBD complex is formed), they contribute with a higher amount of negative charge to the electrolyte than in the idle state (ACE2 receptor only). As a consequence, as α augments one would expect an increase in the output current in the p-branch due to a higher hole concentration induced by the activated receptors and, conversely, a reduction of the output current in the n-branch due to a diminished electron concentration. On the contrary, we observe that IDS decreases (increases) for higher α values in the p-branch (n-branch). A similar behaviour has been also reported in a very recent experimental realization in ref. 7 where the capture of the negatively charged spike protein does not result in a direct increase (decrease) of the p-branch (n-branch) sensor current. In this case the receptor is a spike-protein antigen, positively charged, but the capturing of the negative spike protein should, in principle, reduce the net positive charge at the sensing interface and have the same intuitive consequences above explained.
We can further analyse this behaviour leading our attention to the actual electron and hole densities in the graphene layer for different occupation percentages as a function of VFG (see Fig. 4b). As expected, below VDirac holes are the majority carriers in the channel but their concentration barely changes with α. The inset exhibits a zoom of this p − VFG profiles indicating a slight reduction of the hole concentration as the number of activated receptors increases. On the other hand, above VDirac electrons are dominant, but a non-negligible amount of holes are still present in the n-branch. In addition to this, the amount of electrons in the graphene layer is modulated in a larger extent by α when compared with the holes in the p-branch, so that we can associate both facts with the larger change in IDS observed for the n-type branch. Thus, the electron and hole concentrations, although coherent with the IDS results, still show trends with α that are not intuitive: i.e. the hole (electron) density decreases (increases) as more receptors become activated and therefore more negative charge is present in the electrolyte.
In order to shed light on this issue, we have analysed the longitudinal charge profiles under each receptor with α = 0.6 at two gate biases, one in the p-branch (VFG = −0.5 V) and another in the n-branch (VFG = 0.5 V) (Fig. 5a and b respectively). The charge varies along the channel mimicking locally the changes due to the molecule activation. To analyse in detail these variations, we selected the region under the sensing layer and split it into 10 sub-regions, in correspondence to the 10 receptors, and plot the resulting localized charged profiles, evaluating the impact that a change in the receptor state has on the main carrier distribution (Fig. 5c and d for holes and electrons, respectively). In these plots, dashed (solid) lines correspond to the carrier distribution under activated (idle) receptors. The profiles for each state overlap with small differences associated to the subtleties of the channel local changes. The hole density (Fig. 5c) shows a noticeable increase when the S1RBD (negatively charged) is captured by the ACE2 receptor. Although the peak in the hole concentration under the activated pair is higher, its profile is narrowed and the resulting overall hole concentration (i.e. the integral of these profiles) decreases. For the electron density in the channel (Fig. 5d) the situation is reversed: under the negatively charged receptors the concentration drops significantly (solid lines), but after the complexation of the receptors (dashed lines) the extension of the depleted region becomes narrower, giving rise to a net increase in the electron concentration.
The previous analysis is also supported by the electric field distribution under the molecules. Fig. 5e–h show it for both molecule states (idle at Fig. 5e and f and activated at Fig. 5g and h) and both VFG values in the p- and n-branch (Fig. 5e, g and f, h, respectively). The edges of the graphene layer are indicated by arrow heads aside each plot. The regions where the electric field is normal to the graphene layer correspond to those locations where the carrier concentration changes by a larger extent. When switching from idle (Fig. 5e and f) to activated (Fig. 5g and h) those regions with higher electric field shift from right to left, which is consistent with the charge profiles depicted in Fig. 5c and d. The present analysis evidences the necessity of a fine-grained treatment of the molecular shape and its charge distribution to capture the subtleties of the electrostatic interaction between the receptor and the graphene channel.
Finally, we have determined the absolute and relative change of the GBioFET output current for different values of the receptor occupation (Fig. 6a and b respectively). We opted for this read-out in the sensor response as our results depict high changes in both the p- and n-branch currents. It is worth noting, however, that the contact resistance might saturate and obscure these changes in the on-state current, and the variation in VDirac can alternatively be exploited to define the sensor response in such a case. The sensor presents a noticeable sensitivity and ability to detect SARS-CoV-2 when it is bound to the corresponding ACE2 receptor which assures its specificity. The sensitivity is higher in the n-branch than in the p-branch, reaching relative changes in the output current of 30% in the former case when all the molecules are activated. Its dependence on the gate bias also varies in both cases: while in the p-branch it reaches a maximum (around VFG = 0 V), then decreases and eventually saturates; in the n-branch the sensitivity increases monotonically with VFG indicating that higher positive biases favour the sensing capabilities of the GBioFET. Indeed, the maximum sensitivity within the range of biases simulated here occurs for VFG = 0.7 V. The inset shows the sensor response at this particular bias as a function of the occupation percentage. Interestingly, the sensitivity increases linearly with α, also guaranteeing the linearity of the sensor response at this bias point.
Footnote |
† Electronic supplementary information (ESI) available. See https://doi.org/10.1039/d2na00357k |
This journal is © The Royal Society of Chemistry 2022 |