Graphene BioFET sensors for SARS-CoV-2 detection: a multiscale simulation approach.

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. This journal is © The Royal Society of Chemistry.

Introduction

The world health emergency caused by the Severe Acute Respiratory Syndrome Corona Virus-2 (SARS-CoV-2) has originated an urgent necessity to develop early detection tools to limit and mitigate its uncontrolled spread. Current detection methods for SARS-CoV-2 are mainly based on the extraction of biological samples (e.g. nasopharyngeal or oropharyngeal swabs, blood samples, etc.) for their later analysis with techniques like serology or real-time Reverse Transcription Polymerase Chain Reaction (RT-PCR).[1,2] Alternative methods are being investigated based on imaging procedures, like Computed Tomography (CT),[2,3] X-ray,[4] or DNA amplification as LAMP (loop-mediated isothermal amplification method).[5,6] However, these techniques require a protocol for sample preparation, hours for obtaining a conclusive outcome or sophisticated facilities. Biosensing devices based on field-effect transistors (BioFETs) represent an advanced technology with potential to overcome all these drawbacks.[7]

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.

Methods

The proposed multiscale simulation platform for BioFET treats the two regions that encompass the device, i.e. the graphene and the electrolyte, making use of different descriptions connected through the Poisson equation solved in a two-dimensional cross-section of the structure.

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:where q is the elementary charge, VE is the Fermi level potential (EF = −qVE), n (p) is the electron (hole) density and μn (μp) is the electron (hole) mobility, which includes electric field dependencies and velocity saturation effects.[27] The electron and hole carrier densities at the channel are determined from the combination of the Fermi-Dirac occupation distribution and the density of states of the 2D crystal, which can be extracted from atomistic simulations, semi-empirical models or analytical expressions. Here, we consider a single graphene layer so that, in the range of energies of interest for electronic transport, close to the Dirac point, the band structure is accurately represented by a linear energy–momentum (E–k) relationship.

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:where, c is the i-th ion bulk concentration, cmax, is the maximum allowed concentration, z is the ion valence and VPMF, is the PMF profile. V stands for the 2D potential profile and Vref is the reference potential of the solution. This platform, including the self-consistent solution of the semiconductor and the electrolyte systems, has been already validated obtaining an excellent agreement with experimental results.[31]

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:6vxx[34]) and open state (PDBe ID:6vyb[34]). 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:6m0j[36]).

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: 1o86 [35] and PDBe ID: 6m0j[36] 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 server[41] (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.

Results

We exploit the previously described approach to study the sensitivity of GBioFETs as sensors for the detection of SARS-CoV-2, capturing the finest details in the interaction between the molecules at the sensing interface and the graphene channel. For this purpose, we have considered a graphene monolayer 500 nm-long on top of a SiO2 substrate. Source and drain regions span 50 nm from the edges of the structure and each of them is covered by a 20 nm thick HfO2 layer. In the 400 nm-long remaining region, a 3 nm thick lipid layer that hosts the receptor molecules is considered. Despite the importance of contact resistance in the context of 2DM-based FETs,[42-46] we have consider ohmic contacts focusing on the intrinsic physics of the sensor and how it is impacted by the molecule features. Simulations are carried out with 10 receptors uniformly distributed along the structure. A liquid electrolyte covering the top surface of the sensor is defined by a 1× PBS (Phosphate Buffered Saline) solution. Finally, a reference electrode is immersed on the liquid electrolyte to set a reference voltage. Fig. 3 depicts a schematic of the considered GBioFET.

Fig. 3

Schematic depiction of the GBioFET with a 500 nm-long graphene layer deposited above a SiO2 substrate. 50 nm source and drain regions (black dashed rectangles) are covered by a 20 nm-thick HfO2 layer, while the 400 nm-long channel is covered by a 3 nm-thick lipid membrane that hosts the receptor molecules. To model the substrate of the sample we considered an electrolyte defined by a 1× PBS solution. A reference electrode in direct contact with the solution is also indicated.

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).

Fig. 4

(a) Transfer response of the GBioFET device for different receptor occupation percentages (α). Inset shows the change of the Dirac point when that activation percentage is modified. All the data were obtained for a constant VDS = 0.1 V. (b) Electron (solid) and hole (dashed) concentrations in the graphene layer as a function of the reference electrode bias (VFG) and occupation percentage α. The inset shows a zoom of the hole concentration to highlight the trend of this parameter as α is modified.

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.

Fig. 5

Longitudinal hole (a) and electron (b) profiles in the α = 0.6 scenario. The region under the sensing layer (delimited by red dashed lines) is split into 10 subregions (indicated by dashed black lines) to extract the profiles showed in figures (c) and (d). They correspond to the electron (c) and hole (d) profiles under each receptor that are plotted superimposed to illustrate the changes when the state of the receptor switches from idle (solid lines) to activated (dashed lines). Bottom figures corresponds to the electric field distribution under two receptor molecules in idle state (e, f) and activated state (g, h) for two different gate biases −0.5 V and 0.5 V. They show how the charge of the molecule is redistributed giving rise to a change in the electric field in the graphene layer, the limits of which are indicated by arrow heads. The regions where the vertical component of the electric field is stronger correspond to those where the carrier concentration is modified in a larger extent. When moving from one state to other, the modification of the region with a high vertical electric field gives rise to the changes in the profiles depicted in (c, d).

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.

Fig. 6

Absolute (a) and relative (b) changes of the output current with respect to the 0% occupation scenario. The insets show the response in each case at the bias where the device output shows the largest changes.

Conclusions

We present a pioneering multiscale numerical characterization of graphene-based BioFETs incorporating the description of the charge distribution of the ACE2 receptor and of the ACE2-S1RBD complex that is formed after the virus spike protein anchors to the ACE2 receptor on the sensor surface. The molecular charge distribution is mapped in atomic detail and incorporated into macroscopic device level simulations to assess the GBioFET response to the presence of the target analytes on the sensor surface. The fine-grained level of description proposed here, which is possible due to a multiscale approach to the electrolyte device physics, results in an exceptional level of detail to capture the sensor response. The computational cost at the device level is still high and further improvements would be required to consider larger structures, like the complete spike protein or even the virus. Furthermore, the upgrade of the whole platform towards 3D simulations would be desirable to avoid the projection mapping of the molecular charge into the device configuration. Despite these initial limitations, this work demonstrated the potential of GBioFET to achieve a reliable detection of SARS-CoV-2, particularly when it operates in the n-branch. A detailed study of the electron and hole density variations due to the electric field generated by the ACE2 and the ACE2-S1RBD pair charge distributions, namely when the receptors are idle or activated, enabled us to explain what is, in principle, a counter-intuitive behaviour of the sensor response, that would be obscured by more simplistic computational treatments, rationalizing the observed behaviour in experimentally fabricated sensors. In this case, the ACE2-S1RBD pair has been considered, but future works can explore the combination of this receptor with other molecules to evaluate the impact of these changes in the response of the device. This can be also extended to any other receptor-target pair, and in particular to the investigation of sensor selectivity. Although a selectivity study would require a previous dedicated biochemical investigation of the relevant receptor-target pairs, we believe that the proposed platform constitutes a powerful tool to close the gap between the theoretical predictions and experimental selectivity measurements as it reduces the equivalence between molecules at the device computational level to those with similar electric charge distributions. In summary, we reported a multiscale approach that advances the state-of-the-art of computational tools for the study of biosensors, and which can be exploited in the design of future sensors, including, for example, the analysis of different 2DMs and oxides along with atomic level particularities of the material[47] or the assessment of the impact of the device geometry on its sensitivity. Moreover, the present study might serve as a basis for further detailed investigations, including the impact on the sensor response non-idealities like charge traps,[48] contact resistances, or defects in graphene.

Conflicts of interest

There are no conflicts to declare.