Mario Par

Araújo Luis

Rocha

Dirigée Par

Benoît Piro

Vincent Noël

Isabelle Mme

Navizet

MCU Mme Christelle Yeromonahos

M Florent Barbault

M Vincent Noël

M Benoît Piro

Joaquim E Conceição

Minha Irmã Priscilla

Noivo Fabrice

Eld

Felipe

Fabio Marco

Flavio

Guilherme Dedico

Mario Araujo-Rocha

Florent Barbault

Molecular modelling of interfacial conformational switch of immobilized DNA structure upon target binding

Version Française

Les biocapteurs sont des dispositifs analytiques qui combinent la spécificité d'un composant biologique ou biorécepteur, tel que des protéines, des acides nucléiques, des tissus et même des organites ou des cellules entières, à un composant physique transducteur capable de mesurer et de quantifier l'interaction entre le composant biologique et l'analyte. Une approche très puissante consiste à utiliser des surfaces modifiées par des biorécepteurs dont les propriétés physicochimiques dépendent de l'environnement voisin.

Dans ce travail, nous effectuons des études computationnelles sur un biorécepteur à base d'ADN pour un biocapteur de type EGOFET en utilisant différentes méthodes de calcul. Nous avons choisi l'ampicilline (AMP) comme molécule d'intérêt pour le biocapteur. Afin de détecter la présence d'AMP avec notre biorécepteur, nous avons ajouté une séquence d'aptamère anti-AMP à notre séquence globale d'ADN, dans le but de mieux comprendre le processus de transduction en analysant les changements de distribution ionique dans notre système sans ou avec l'AMP. Nous montrons nos résultats sur la construction de notre structure de bioreconnaissance modèle à la fois dans les états off et on, représentant le système en l'absence et en présence d'ampicilline (AMP), respectivement.

Nous avons effectivement trouvé des conformations stables pour les deux états, y compris l'élucidation de la structure complexe AMP-aptamère par docking moléculaire. Nous concluons que notre système proposé est prometteur pour tout type de transduction électrique basée sur un changement de capacité ou d'effet de champ, et nous avons élucidé certains aspects du processus de bioreconnaissance qui ne seraient pas si faciles à réaliser expérimentalement, tels que les profils ioniques et les mécanismes de déshybridation possibles. Mots clés : Biocapteurs ; Transistor à effet de champ ; Dynamique Moléculaire ; ADN ; Aptamères ; Docking Moléculaire.

Abstract -English Version

Biosensors are analytical devices that combine the specificity of a biological component or bioreceptor, such as proteins, nucleic acids, tissues and even entire organelles or cells, to a physical transducing component capable of measuring and quantifying the interaction between the biological component and the analyte. One very powerful approach is to use bioreceptor-modified surfaces whose physiochemical properties depend on the neighboring environment.

In this work, we perform computational studies on a DNA-based bioreceptor for an EGOFET-type biosensor using different computational methods. We have chosen ampicillin (AMP) as our molecule of interest for the biosensor. In order to detect the presence of AMP with our bioreceptor, we have added an anti-AMP aptamer sequence to our global DNA sequence, with the goal of better understanding the transduction process by analyzing the changes in ionic distribution in our system in both absence and presence of AMP.

We have shown our results on the construction of our model biorecognition structure in both off-and on-states, representing the system in the absence and presence of ampicillin (AMP), respectively. We have effectively found stable conformations for both states, including the elucidation of the complex AMP-aptamer structure. We conclude that our proposed system is promising for any kind of electrical transduction based on a change in capacitance or field effect, and we have elucidated certain aspects of the biorecognition process that would not be so easily done experimentally, such as the ionic profiles and the possible de-hybridization mechanisms. Keywords: Biosensors; EGOFETs; Molecular Dynamics; DNA; Aptamers; Molecular Docking. I would like to thank my supervisors who guided and supported me until the very end of this project: Professor Benoît Piro, I could not have asked for a more patient, understanding and dedicated supervisor. Professor Vincent Noël, who never refused a good scientific discussion and his passion for science is a true inspiration for me. Professor Florent Barbault, he is the one that instigated my passion for Molecular Dynamics. It is not an exaggeration to say this thesis would not have existed without him. I would also like to thank all the Professors from the CTM team: Mahamadou Seydou, Eric Bremmond, François Maurel and Antonio Monari. The CTM team is an oasis of healthy and pleasant work environment thanks to them. I'd also like to thank my fellow PhD students in the team Seydou Sall and Sylvain Pitié (who was also the best karaoke partner I have ever had), as well as Roland Nguyen, who gave me TI support uncountable times during the thesis. I would also like to thank the entire BIOSS team, Professors Guillaume Anquetin, Nicolas Battaglini, Samia Zrig and Steeve Reisberg. Fellow PhD students Justin Lemarchand, Nga Dau, Baptiste Rouchié, Rassen Boukraa, Samia Mekhmoukhen, as well as former PhD student Sandra Vasilijevic and the newcomers Clara dos Santos and Pierre Guermonprez who unfortunately I didn't get to spend enough time with. I would also like to thank former members Khalil Chennit, Natalie Bridonneau, Justine Quinet and Flavien Mouillard. A special thanks to Professor Giorgio Mattana for the great theoretical discussion on the electrically charged surfaces and for indicating the approximate location of Professor Vincent Noël on a regular basis.

At ITODYS, I have met great people who made my days there very pleasant. I would like to thank Larissa Chaperman, Sam Hoff, Christine Ranjan, Célia Sahli, Sandos Darwish, Hedi Kouki, Imen Hnid, Asma Torkhani, Chaima Tounakti, Yahya Hami, Mathieu Bastide, Yara Al Chami, Mingyang Liu, Xinlei Yao and Xiaonan Sun. I would also like to thank the people in the administration of both the lab and the UFR Chimie, with a special thanks to Pierre-François Quenin, Juliette Le, Juliette Leblanc and Philippe Cabaniols.

List of Abbreviations

General introduction

Biosensors are analytical devices that combine the specificity of a biological component or bioreceptor, such as proteins, nucleic acids, tissues and even entire organelles or cells, to a physical transducing component capable of measuring and quantifying the interaction between the biological component and the analyte. For a simple utilization of a biosensor, including by untrained average people, it is preferable for the detection to be done without adding additional reagents nor modifying the original sample (label-free and reagent-less detection). One very powerful approach is to use bioreceptormodified surfaces whose physiochemical properties depend on the neighboring environment. In such conditions, bioreceptor/target complex formation leads to a direct surfaces properties change.

Among techniques that are capable of detecting changes on the environment of modified surfaces, there are the electrolyte-gated organic field effect transistors, or EGOFETs, which have been effectively developed as transducers in literature, being able to directly monitor biorecognition processes between an analyte and its bioreceptor at one of its interfaces of its electrolyte.

In this work, we perform computational studies on a DNA-based bioreceptor for an EGOFET-type biosensor using different computational methods. We have chosen ampicillin as our molecule of interest for the biosensor. Ampicillin (AMP) is a β-lactam antibiotic that has been extensively used for treating bacterial infection in both humans and animals since the 1960s 1 thanks to its activity against gram-negative bacteria, something that penicillin treatments were not capable until then 2 . This extensive usage has caused ampicillin to become a relevant environmental pollutant over time. In order to detect the presence of AMP with our bioreceptor, we have added an anti-AMP aptamer sequence to our global DNA sequence, with the goal of better understanding the transduction process by analyzing the changes in ionic distribution in our system in both absence and presence of AMP.

Biosensors

Biosensors are analytical devices that use the specificity of a biological structure to recognize a molecule of interest. Because biosensors in general have been extensively described and reviewed, we will not give here more details and we invite the reader who wants to know more about biosensors to read some excellent reviews available in the recent literature.

Among other parts, biosensors are composed of a biorecognition module that is connected to a physio-chemical transducing component. The biorecognition module can be a protein, an antibody, a nucleic acid and even entire cells or tissues, while the transducing component is in charge of detecting the interaction between the biorecognition module and the analyte and quantifying it, providing readability for the user (Figure 1). This can be done in many ways, such as optical, piezoelectric, electrochemical or gravimetric methods, for example. The concept of biosensor was first introduced in the seminal work of Clark and Lyons in 1962 describing a glucose biosensor using an oxygen electrode treated with a thin layer of glucose oxidase between two semipermeable membranes 3 . The concept of a glucose biosensor using glucose oxidase as bioreceptor was further developed by Updike and Hicks in 1967, who reported the immobilization and stabilization of the enzyme in a polyacrylamide gel over the electrode 4,5 , until finally the first commercial glucose biosensor was launched in 1975 6 . This biosensor, however, was not portable and its high upfront cost caused it to be exclusively used in clinical applications 6,7 . Nowadays, analytical methods are omnipresent in medicine, however, they usually employ amplification methods that demand a laboratorial treatment by qualified personnel, with examples of such methods being Enzyme-Linked ImmunoSorbent Assay (ELISA) and Real-time Reverse Transcription Polymerase Chain Reaction (rRT-PCR), with such analyses among the main methods of viral detection during the SARS-CoV-2 (COVID-19) pandemic 8 since 2020. However, the use of biosensors (i.e. autonomous and portable devices) in medical laboratories is not common yet. But their use by the untrained population is increasing steadily, outside laboratories, with the example of methods such as the Lateral Flow (LF) that is used in commercial devices to detect pregnancy, HIV and SARS-CoV-2 antibodies 8 , among other analytes 9 .

The glucose biosensor is also nowadays widely available as a portable device for the monitoring of the blood sugar levels for both clinical usage and home care, and it is the pioneer on what today is defined as point-of-care testing: medical diagnostic testing that can be done directly at the patient's location instead of sending the biological sample to an external laboratory. Point-of-care testing has been a staple in medicine 10 , reducing general hospital stay time and improving treatment adherence 11 . A study in an intensive care unit in 2004 showed that point-of-care testing had decreased mortality in neonates after congenital heart disease surgery between 2001 and 2003 from 12% to 4.3% 12 .

In the recent decades, the integration of point-of-care testing to every day's objects has been increasingly discussed. The idea of coupling biosensors to smartphones has gained interest in literature [13][14][15][16] . Smartphones have great potential as biosensing devices as they usually contain a screen for display and user interface, some detectors of external signal, such as cameras 17 , light sensors 18,19 and even the microphones 20 that can be used to record the biosensor signal as a transducer, multiple communication pathways (Internet, Wi-Fi, Bluetooth, NFC) as well as being potential energy sources.

One other increasingly popular way of integrating biosensors to every day's life is in the form of wearable sensors 21,22 . For example, the use of programmed bacteria in wearable patches to detect certain chemical compounds on the skin has been reported by Liu et al 23,24 between 2017 and 2018. Wang et al 25 similarly published a work in 2017 on a bacterial bio-hybrid material used to detect moisture on skin coupled with smart flaps that open when the moisture is high to create smart clothing that aerates with the heating up of the body. Wearable glucose meters are also constantly in development 26 aimed at the measurement in different body fluids such as urine 27,28 , interstitial fluids [29][30][31][32] , saliva 33,34 , tears 35,36 , serum 37 or sweat 38,39 .

Other important aspects for the easiness of usage by the untrained general population are the notions of reagent-less and label-free detection. Reagent-less detection refers to a biosensor capable of functioning by the insertion of the material to be analyzed without the further addition of any other chemical reagent, while label-free detection refers to the analysis of the material to be analyzed without any previous chemical modification of this material. These two concepts together eliminate two steps of the detection that can be a source of errors when performed by the untrained population.

Coupling biosensors to smartphones and wearables, however, often requires a certain degree of miniaturization of the biosensor to portable sizes, to fit in systems that are called "Lab-on-a-Chip", a small device incorporating the processes involved in a laboratory analysis over an integrated circuit unit. These devices work with very small fluid volumes and, to allow the detection in suboptimal body fluids such as interstitial fluids, tears or saliva, there is a demand for lower limits of detection in biosensors, reaching in certain cases the necessity of single molecule detection.

Single molecule detection is a concept that became increasingly popular since the 90s 40 , with seminal works such as the first imaging of single fluorophores on a surface 41 and in solution 42 . In 2007, label-free single molecule detection of interleukin-2 in serum was reported by Armani et al 43 in optical silica microcavities of diameter 80 μm by surface plasmon resonance, with concentrations down to 300 aM being detected.

In 2018, Macchia et al 44 reported a label-free single molecule detection using a millimeter-sized electrolyte-gated organic field effect transistor (EGOFET). In our study, we were interested in the EGOFET configuration for our proposed biorecognition structure, and we will further discuss this type of device in the next section of this chapter.

The COVID 19 pandemic also pushed further studies in effective detection of the SARS-CoV-2 biomarkers. The detection of SARS-CoV-2 biomarkers in healthy individuals is particularly critical due to the lack of symptoms presented in a large part of the infected population, this population being able to contaminate further individuals. In March 2021, Meena et al 45 reported the multiplexed detection of both antigen and RNA of SARS-CoV-2 from nasopharyngeal swab samples to a single molecule level by fluorescence with the help of an optofluidic chip.

In July 2021, Guo et al 46 reported the development of a device based on an organic electrochemical transistor capable of single molecule detection of SARS-CoV-2 and MERS antigens in saliva, illustrated on Figure 2. Organic electrochemical transistors are structurally very similar to EGOFETs, with the major difference being the ions present in the electrolyte being able to penetrate into the channel instead of accumulating near the surface as in EGOFETs. As bioreceptor, they have used a Nanobody® 47 agent sensitive to a spike protein present at the virus capsid. The authors believe that their system could be used for early detection of viral load in saliva of asymptomatic individuals, which are the main spreaders of the virus 48,49 . But what is the physical limit of single molecule detection (and its miniaturization)?

In the seminal theoretical work of Sheehan and Whitman in 2005 50 , the authors argue that a sub-femtomolar detections in nanometric detectors should take from hours to days to take place if we take into account the diffusion of the molecules in solution, even though the previously mentioned work by Armani et al 43 indicated otherwise. In 2008, intrigued by this discrepancy between theory and experimental results, Squires, Messinger and Manalis 51 published a very thorough theoretical study on the optimization of a biosensor's dimensions and their insertion on microfluidic systems to come around the previous study results of Sheehan and Whitman. But even then, the authors admitted that there were systems that did not follow their "speed limit" for mass transport diffusion even taking into account electrostatic effects, especially those involving nucleic acids [52][53][54] .

If a large part of the future of biosensing devices is in miniaturization and single molecule detection, the physical limits of detection eventually will have to be dealt with. This could be done by maintaining a macroscale surface area of detection or by improving association rates between analyte and bioreceptor. However, biosensors can only be miniaturized to a certain point without also decreasing the area of detection, and at this point, improving the association rates becomes essential.

It is, therefore, highly necessary to better understand the molecular interactions that may improve the detection, as well as the mechanisms behind the transduction between the biorecognition processes and the physio-chemical signals to be measured. This work focuses on the latter, in an attempt to unveil the transduction process of an EGOFET-type biosensor, which is the topic of the following section.

EGOFETs -Theory and applications in biosensors

The term EGOFET stands for Electrolyte-Gated Organic Field-Effect Transistor, and it describes an organic field-effect transistor in which the dielectric component is an electrolyte. An EGOFET is composed of three terminals (source, drain and gate), an organic semiconductor (OSC) and the electrolyte. When a voltage is applied at the drain terminal (𝑉 𝐷 ), the resulting current (𝐼 𝐷 ) can be modulated by applying a tension in the gate terminal ( 𝑉 𝐺 ). This is achieved by a field effect between the gate and the organic semiconductor, through the electrolyte, which causes the accumulation or depletion of the charge carriers in the semiconductor, therefore increasing or decreasing its conductivity.

This structure is further illustrated on Figure 3. For a p-type semiconductor (examples: PEDOT:PSS, DPP-DTT) 55 , a negative 𝑉 𝐺 causes an increase of charge carriers near the surface of the semiconductor while a positive one depletes the surface of charge carriers, while for a system containing a n-type semiconductor (examples: naphthalene tetracarboxylic diimides, azulene derivatives) 55 the opposite is true. However, in both cases the applied 𝑉 𝐺 on the system needs to surpass in absolute value a certain threshold in order to properly turn on the transistor; this minimal bias is called threshold voltage (𝑉 𝑇 ).

The gating (appliance of voltage on the gate electrode terminal) of EGOFETs leads to the formation of two electrical double layers, one at the interface between the gate and the electrolyte and another at the interface between the electrolyte and the semiconductor.

Both layers have high capacitance (even if the gate/electrolyte interface is often higher); therefore it is possible to run EGOFET systems using low voltages, with less than 1 V 56 .

As the resulting currents are also dependent on the capacitance of those two double layers, EGOFETs are intrinsically sensitive to changes on the electric environment in these layers. This can be exploited to create EGOFET-based biosensors via functionalization of one of the transistor interfaces (gate or channel). Any change of the ionic distribution in the doubles layers sufficient to impact their capacitance can be amplified thanks to the field effect, i.e. a small capacitance change leads to a large output current variation. Some results, however, require a further analysis of the electrical environment near the bioreceptors in EGOFETs. In 2015, Palazzo et al 57 reported the detection of interactions beyond the Debye's length using EGOFETs, meaning that events occurring outside the double layer can be also monitored. The Debye length is a common measure of a charge's effect in a solution and how far its electrostatic effects persist until being screened out by a magnitude of 1/e, and it is also commonly considered a classical limit for electrochemical screening, that showed to be problematic in other FET-type systems 58,59 . The Debye length (𝜆 𝐷 ) in electrolytes is usually defined by the following equation:

𝜆 𝐷 = √ 𝜀 0 𝜀 𝑟 𝑘 𝐵 𝑇 ∑ 𝑞 𝑖 2 𝑐 𝑖 0 𝑖 2 (1)
being 𝜀 0 the dielectric permittivity in vacuum, 𝜀 𝑟 being the relative permittivity of the solvent, 𝑘 𝐵 the Boltzmann constant, 𝑇 the temperature in Kelvin, 𝑞 𝑖 and 𝑐 𝑖 0 being the charge and concentration of the ion 𝑖. One can attempt to explain the results by Palazzo et al 57 The behavior of an EGOFET is usually characterized using two types of graph curves.

The first one, called output curve indicates the behavior of 𝐼 𝐷 with respect to 𝑉 𝐷 under different 𝑉 𝐺 values, while the transfer curve shows the behavior of √𝐼 𝐷 under a fixed bias on a saturated 𝑉 𝐷 with respect to the VG. In this work, we will only discuss the general details of the output curves. Figure 5 shows a generic example of an output curve (obtained on devices fabrication in the group). Normally, the output curve follows two types of regime: linear regime and saturated regime. The linear regime occurs when the applied 𝑉 𝐷 ≪ (𝑉 𝐺 -𝑉 𝑇 ) , and it is ohmic in nature, with respect to 𝑉 𝐷 . The equation that describes the linear regime is the following, with 𝑊 being the channel's width, 𝐿 its length, 𝐶 𝑖 being the total dielectric capacitance per unit area and μ the field-effect mobility of the charge carriers:

𝐼 𝐷 = 𝑊 𝐿 𝐶 𝑖 μ (𝑉 𝐺 -𝑉 𝑇 )𝑉 𝐷 (2)
The saturated regime happens when 𝑉 𝐷 becomes so strong that the charge carriers in the semiconductor channel start to be depleted at the drain contact. This causes a stabilization of the current along the channel, reaching a saturation plateau with 𝐼 𝐷 following the following equation:

𝐼 𝐷 𝑆𝑎𝑡𝑢𝑟𝑎𝑡𝑒𝑑 = 𝑊 2𝐿 𝐶 𝑖 μ (𝑉 𝐺 -𝑉 𝑇 ) 2 (3)
One more important detail to remember is that 𝐶 𝑖 is dependant of the capacitances of the double layers in the gate/electrolyte and electrolyte/semiconductor interfaces by the following equation:

𝐶 𝑖 -1 = 𝐶 𝑂𝑆𝐶/𝐸𝑙𝑒𝑐 -1 + 𝐶 𝐺𝑎𝑡𝑒/𝐸𝑙𝑒𝑐 -1 (4)
In both regimes, the final 𝐼 𝐷 value has a direct relation to the total capacitance of the system 𝐶 𝑖 , and this is the relation to be exploited when functionalizing the interfaces with bioreceptors.

EGOFETs have been vastly studied in the last decade: to this day, a total of 392 articles are listed on Google Scholar citing EGOFETs from 2011 to 2021. The articles reported in this section are only a small fraction of the available literature produced in the recent years.

Our team is deeply involved in the development of EGOFETs. Among contributions to the academic discussion on the topic, in 2010, Kergoat et al [START_REF] Kergoat | A Water-Gate Organic Field-Effect Transistor[END_REF] published a work on the development of the first (technically not EG)OFET device gated with pure water using a single crystal of rubrene as semiconductor. In 2017, the use of EGOFETs as an immunoassay for bisphenol A detection was reported by Piro et al, by co-crystallizing an alkyl derivative of the analyte with the organic semiconductor 61 . By adding a bisphenol-Asensitive antibody, any bisphenol A in the electrolyte would cause a competitive dissociation of the antibody in the semiconductor's interface, therefore changing the capacitance of the double layer and being detected as a change in 𝐼 𝐷 . The authors explain that this method could be extended to the detection of many different small molecules. In 2018, the same group 62 reported a pH-sensitive EGOFET by the functionalization of a gold gate electrode with a switchable hydrogel of poly(acrylic acid). This hydrogel could increase or decrease its volume (and its thickness) according to the pH in the electrolyte, again changing the capacitance of the gate interface system and the resulting 𝐼 𝐷 .

To reach the single molecule detection, we should master all chemical events at the atomic scale and, hence, this deciphering can be undertaken with molecular modelling.

As quickly discussed in the previous section, single molecular detection of immunoglobulin-G in EGOFETs was reported in 2018 by Macchia et al 44 On December 2020, Sensi et al 63 published an EGOFET-based biosensor integrating a biological drug (Nivolumab, a human monoclonal antibody indicated for immunotherapy of certain types of cancer) to the gate interface to detect anti-drug antibodies, which signalizes an immune response against the drug, a main issue in immunotherapy effectiveness 64 . The authors report a limit of detection of 100 fM for this biosensor.

On March 2021, the detection of miRNA-21 at picomolar levels was reported by Selvaraj et al [START_REF] Selvaraj | Label Free Detection of MiRNA-21 with Electrolyte Gated Organic Field Effect Transistors (EGOFETs)[END_REF] using a double-gated EGOFET with one gate functionalized with 2mercaptoethanol and the other functionalized with a thiolated RNA probe that is complementary to the target sequence. The authors also estimated a binding free energy 20% lower for the hybridization of the strands grafted to the gate electrode compared to solution, which they assumed to be either from an inhibiting effect by the surface or competing reactions. This effect on binding free energy next to the interface is of great relevance to our studies.

While technically not an EGOFET, it is important to mention the usage of an electrolyte-gated graphene field effect transistor-based biosensor for the detection of SARS-CoV-2 in nasopharyngeal swab fluids reported by Seo et al 66 in April 2020 with limit of detection of 242 copies per mL in clinical samples, by the functionalization of the graphene channel with anti-spike-protein antibodies.

Biosensors using aptamers as bioreceptors are also called aptasensors. The notion of aptamer is further discussed in this chapter, as well as some general aspects of nucleic acids in general in section 4. To date, the only effective use of aptamer nucleic acid sequences in EGOFETs to detect its molecule of interest was reported by Pallu et al [START_REF] Pallu | A DNA Hydrogel Gated Organic Field Effect Transistor[END_REF] , using a DNA hydrogel containing a sequence sensitive to vascular endothelial growth factor VEGF (a cancer biomarker [START_REF] Qureshi | Capacitive Aptamer-Antibody Based Sandwich Assay for the Detection of VEGF Cancer Biomarker in Serum[END_REF] ). This was done by using the bulk DNA hydrogel as a solid electrolyte gating the transistor. Outside this usage, the only other published article describing an EGOFET aptasensor is a work by Berto et al [START_REF] Berto | EGOFET Peptide Aptasensor for Label-Free Detection of Inflammatory Cytokines in Complex Fluids[END_REF]

Nucleic acids

Nucleic acids are a class of biopolymers found in every type of known life form.

There are two types of nucleic acids found in nature: deoxyribonucleic acid, more commonly known as DNA, and ribonucleic acid or RNA. While other types of synthetic nucleic acids are also commercially available, in this section we will only discuss the natural occurring types of nucleic acids, with a particular focus on DNA.

DNA and RNA are involved in a series of biological processes that involve information transmission, such as storage and processing of genetic information, transduction and synthesis of proteins. This is thanks to their capacity of codification through the sequencing of their monomers, also called nucleotides. The basic structure of a nucleotide is shown on Figure 7 (left). The right section of Figure 7 shows the structural formula of all these nucleobases.

One important aspect of these nucleobases is that they have H-bond acceptor and donor sites. These are the main interaction sites between two bases involved in the canonical base pairing that we find in both DNA (such as in double helixes) and RNA (hairpins, intra-strand double helixes). This pairing generally follows the rules of Watson-Crick 70 due to the number and position of these sites (see Figure 8): adenine pairs well with thymine (or uracil) as both have one acceptor and one donor H-bond site at similar distances from each other but at inversed positions, therefore generating an optimal pairing where one base's donor site interacts with the other's acceptor site. The same is true for the cytosine-guanine pairing, whereas with three H-bond sites: guanine has the motif acceptor-donor-donor in its H-bond sites while cytosine has the opposite motif, donor-acceptor-acceptor. This alternation of donor and acceptor sites also induce a favorable long-range electrostatic interaction. However, non-Watson-Crick base pairs are also possible, although generically less energetically favorable. Methods of propagation such as Polymerase Chain Reaction (PCR) 72 have opened up the possibilities of large-scale DNA synthesis in laboratory, while the synthesis by the nucleoside phosphoramidite method 73 permits the chemical de novo synthesis of oligonucleotides of up to 200 bases and lately there were reports of highly efficient methods of iterative nucleobase addition [START_REF] Palluk | De Novo DNA Synthesis Using Polymerase-Nucleotide Conjugates[END_REF] . This sort of manipulation of DNA sequences opens up other possible utilities for the nucleic acids outside their biological use, although in this case, DNA is usually preferred over RNA due to its higher stability.

In our work, we are interested in another non-biological aspect of the versatility of nucleic acids: the development of sequences with a high affinity to an analyte. These oligonucleotides are called aptamers, and they are possible due to the high flexibility of nucleic acids with a high modularity of available H-bond donor and acceptor sites, allowing the oligonucleotide to adopt a conformation on which the interaction is maximized.

Aptamers and aptasensors

Although the word aptamer can express either a nucleic acid or an oligopeptide that is designed for high binding affinity to a target, in this section we will only discuss aptamers based on nucleic acids. Since their first appearances in the 90s in scientific literature and the seminal work of Drolet et al [START_REF] Drolet | An Enzyme-Linked Oligonucleotide Assay[END_REF] in 1996 using a RNA aptamer for the detection of human vascular endothelial growth factor, aptamers have become a topic of extreme relevance in domain of biosensors and clinical diagnostic [START_REF] Stanciu | Recent Advances in Aptamer-Based Biosensors for Global Health Applications[END_REF] .

Aptamer sequences can be experimentally obtained by a method called SELEX (Systematic Evolution of Ligands by EXponential enrichment) [START_REF] Tuerk | Systematic Evolution of Ligands by Exponential Enrichment: RNA Ligands to Bacteriophage T4 DNA Polymerase[END_REF] . This method consists of an initial synthesis of a random pool of oligonucleotide sequences of up to 10 15 different sequences [START_REF] Ellington | In Vitro Selection of RNA Molecules That Bind Specific Ligands[END_REF] , followed by a selection process in which the random pool is incubated with immobilized target molecules. The best binding sequences are then extracted and amplified by a PCR step that allows for random mutations on these sequences, creating the next generation pool of oligonucleotides for the incubational step. The iteration of many of these steps leads to the selection of an oligonucleotide that has a high binding affinity to the target molecule (Figure 9). Aptamers have been described as highly sensitive to many different types of target proteins, such as thrombin [START_REF] Deng | Aptamer Binding Assays for Proteins: The Thrombin Example-A Review[END_REF] and vascular endothelial growth factor (VEGF) [START_REF] Kaur | Probing High Affinity Sequences of DNA Aptamer against VEGF165[END_REF] , as well as targeting small molecules, such as cocaine [START_REF] Stojanovic | Aptamer-Based Folding Fluorescent Sensor for Cocaine[END_REF] , tetracycline [START_REF] Förster | Conformational Dynamics of the Tetracycline-Binding Aptamer[END_REF] , adenosine/ATP [START_REF] Huizenga | A DNA Aptamer That Binds Adenosine and ATP[END_REF] , and 8-oxo-dG [START_REF] Choodet | Combined in Silico and in Vitro Study of an Aptasensor Based on Citrate-Capped AuNPs for Naked-Eye Detection of a Critical Biomarker of Oxidative Stress[END_REF] , a modified deoxyguanosine that is biomarker of oxidative stress.

Aptamers have also been shown to be capable of single molecule detection when coupled to nanostructures such as nanopores [START_REF] Lin | Selective Single Molecule Nanopore Sensing of Proteins Using DNA Aptamer-Functionalised Gold Nanoparticles[END_REF][START_REF] Sze | Single Molecule Multiplexed Nanopore Protein Screening in Human Serum Using Aptamer Modified DNA Carriers[END_REF] or single molecular junctions with carbon nanotubes [START_REF] Liu | Single-Molecule Detection of Proteins Using Aptamer-Functionalized Molecular Electronic Devices[END_REF] and graphene [START_REF] Chen | Ultrasensitive Detection and Binding Mechanism of Cocaine in an Aptamer-Based Single-Molecule Device[END_REF] , using fluorescent colocalization methods [START_REF] Zhang | Single Molecule Fluorescent Colocalization of Split Aptamers for Ultrasensitive Detection of Biomolecules[END_REF] , or using a single molecule bioreceptor [START_REF] Cervantes-Salguero | Single-Molecule DNA Origami Aptasensors for Real-Time Biomarker Detection[END_REF] .

As previously mentioned in the section 3, the only published article known to date reporting an EGOFET nucleic aptasensor is the one by Pallu et al [START_REF] Pallu | A DNA Hydrogel Gated Organic Field Effect Transistor[END_REF] , published in 2019.

First, the authors report the usage of a DNA hydrogel directly in contact with the channel as gating electrolyte, however the output current signal is of lesser quality than adding the hydrogel to a drop of Phosphate-Buffered Saline (PBS) over the channel, with less field effect, high hysteresis, and low stability due to the drying of the hydrogel (Figure 10 a,b).

Using the configuration with the PBS drop, the authors reported the detection of vascular endothelial growth factor VEGF, by the means of following the evolution of the output drain current over time, which increases upon recognition of the aptamer (even if the response time is long due to diffusion limitations) and, comparatively, stays constant with a nonfunctionalized hydrogel (Figure 10, bottom). In the end, the authors estimated a total of 175 ± 22 pg of the target molecule bound to around 50 mg of DNA hydrogel. In this thesis, we are interested in an aptamer sequence that was first described by Song and collaborators in 2012 [START_REF] Song | Aptasensor for Ampicillin Using Gold Nanoparticle Based Dual Fluorescence-Colorimetric Methods[END_REF] . It is specific to ampicillin (AMP), a beta-lactam antibiotic.

In fact, this paper describes three DNA aptamer sequences that are sensitive to AMP found by the SELEX method: AMP4 (5′-CACGGCATGGTGGGCGTCGTG-3′), AMP17 (5′-GCGGGCGGTTGTATAGCGG-3′) and AMP18 (5′-TTAGTTGGGGTTCAGTTGG-3′, the sequence we have chosen for our work), all considered with good sensitivity and specificity by the authors, with estimated dissociation constants of 9.4 nM, 13.4 nM and 9.8 nM, respectively [START_REF] Song | Aptasensor for Ampicillin Using Gold Nanoparticle Based Dual Fluorescence-Colorimetric Methods[END_REF] , and low interaction even with similar beta-lactam antibiotics, such as amoxicillin and benzylpenicillin. The predicted secondary structure of these three sequences is shown on Figure 11. The authors have theorized that the common regions -GGT(T)found in the loop with a CG base pairing at the joint of the loop with the stem played a role in the binding of AMP. In the same article, the effective usage of the sequence AMP17 was reported for dual fluorescence-colorimetric aptasensor using gold nanoparticles [START_REF] Song | Aptasensor for Ampicillin Using Gold Nanoparticle Based Dual Fluorescence-Colorimetric Methods[END_REF] . The complete essays using the other two sequences were not reported by the authors, as an initial test showed that AMP17 was showing slightly more sensitive results for this type of sensor.

In 2019, another AMP aptasensor, this time using the AMP17 sequence, was reported by Taghdisi et al [START_REF] Taghdisi | An Electrochemical Sensing Platform Based on Ladder-Shaped DNA Structure and Label-Free Aptamer for Ultrasensitive Detection of Ampicillin[END_REF] . This sensor is built upon a "ladder-shaped" DNA structure that is supported by the aptamer, as seen on Figure 12, extracted from the article. The authors report a detection limit of 1 pM with a linear range from 7 pM to 100 nM and high selectivity towards other similar antibiotics. There are a few things to improve in their system, however. Their sensor depends on a redox reagent in the form of Fe(CN)6 4-/3-and the DNA structure is not regenerative after the detection. On top of that, there is no amplification of signal: the detection is done by a single electron interaction, which is not a problem in itself, but can be improved using an EGOFET configuration. Overall, their results are very promising and one very relevant point to us from this work is the ability of the aptamer AMP17 to go from a hybridized state to a folded recognition hairpin in the presence of AMP. Even though we do not use the same sequence in this work, we expect the sequence AMP18 to be able to do the same.

The choice of the aptamer sequence in this thesis came instead from a published article by Yu and Lai [START_REF] Yu | A Reagentless and Reusable Electrochemical Aptamer-Based Sensor for Rapid Detection of Ampicillin in Complex Samples[END_REF] in 2018 in which they used AMP18 in a reusable, electrochemical affinity aptasensor. In their device, the aptamer sequence was grafted to a gold electrode treated with 6-mercapto-1-hexanol (MCH) from the 5' terminus while containing a methylene blue probe at 3'. The transduction would come from a conformational switch at the aptamer in the presence of AMP and this switch changes the average distance between the probe and the electrode, which can be detected (Figure 13). Two important things to extract from their result is the effectiveness of the aptamer sequence in binding to the AMP molecule even when grafted to a surface with only a sixcarbon spacer and how the AMP molecule induces the conformational switch in the aptamer. Aptamers in general are known for losing their binding affinity when they are too close to the grafting point on surfaces [START_REF] Balamurugan | Surface Immobilization Methods for Aptamer Diagnostic Applications[END_REF] , possibly by the loss of chain mobility rendering the effective binding more difficult, so this sequence's ability to maintain a high affinity to AMP when grafted this close to the surface is very interesting. The authors tested the responses of the sensor to similar antibiotics, with results indicating high selectivity towards AMP.

To conclude this section, even though aptamers have been successfully used in biosensors with high sensitivity and selectivity, there is only a single report of the usage of nucleic acid aptamers in EGOFETs in literature [START_REF] Pallu | A DNA Hydrogel Gated Organic Field Effect Transistor[END_REF] . This was done by using a very large amount of aptamer-functionalized DNA hydrogel (50 mg) even though there are many reports of sensitive aptasensors using smaller DNA structures [START_REF] Cervantes-Salguero | Single-Molecule DNA Origami Aptasensors for Real-Time Biomarker Detection[END_REF][START_REF] Song | Aptasensor for Ampicillin Using Gold Nanoparticle Based Dual Fluorescence-Colorimetric Methods[END_REF][START_REF] Taghdisi | An Electrochemical Sensing Platform Based on Ladder-Shaped DNA Structure and Label-Free Aptamer for Ultrasensitive Detection of Ampicillin[END_REF][START_REF] Yu | A Reagentless and Reusable Electrochemical Aptamer-Based Sensor for Rapid Detection of Ampicillin in Complex Samples[END_REF] . This work aims to clarify the transduction process of an EGOFET aptasensor in hopes of establishing theoretical groundwork for more efficient future EGOFET aptasensors.

Ampicillin as a pollutant

Ampicillin is a beta-lactam antibiotic first discovered in 1958 [START_REF] Iupac; Fischer | Analogue-Based Drug Discovery[END_REF] that has been extensively used since the first description of its pharmacology in the 1961 by Acred et al 1 and its introduction to market in 1964 96 , thanks to its effect on gram-negative bacteria, on which penicillin treatments were not effective at the time 96 . Its structural formula is represented in Figure 14. Ampicillin is effective on many different types of infections, such as bacterial meningitis, salmonellosis, and general infections by E. coli or bacteria of the genus Staphylococcus 97 . It is considered an essential medicine by the World Health Organization [START_REF]World Health Organization Model List of Essential Medicines: 21st List[END_REF] and used both in veterinary feed [START_REF] Centner | Recent Government Regulations in the United States Seek to Ensure the Effectiveness of Antibiotics by Limiting Their Agricultural Use[END_REF] and humans. However, bacterial resistance to ampicillin has been reported since the 70s in literature 100 and has been a frequent marker of environmental antibiotic resistance [101][102][103][104] . As resistant bacteria pose a serious threat to healthcare 105 , this has motivated governmental agencies to regulate the use [START_REF] Centner | Recent Government Regulations in the United States Seek to Ensure the Effectiveness of Antibiotics by Limiting Their Agricultural Use[END_REF] of antibiotics and therefore, there is a need for more careful antibiotic analyses in general 106,107 . Therefore, the main motivation of choosing ampicillin as the biosensor target in this thesis comes from its environmental importance and the need of sensitive analyses, as well as the existence of well-documented aptamer sequences [START_REF] Song | Aptasensor for Ampicillin Using Gold Nanoparticle Based Dual Fluorescence-Colorimetric Methods[END_REF] that are sensitive to it, as discussed in section 4.1.

State of the art in the atomic modelling of bioreceptors

This is a thesis in computational chemistry in which molecular modelling is applied to a biosensor context. Therefore, it is important to document the most recent advances in the modelling of biosensors, which is the aim of this section, with a focus on molecular mechanics and molecular dynamics studies. These methods will be depicted in Chapter 2.

Computational studies on the interactions bioreceptor-ligand in solution are plentiful in literature [108][109][110][111][112][113][114][115] and there are many methods that are optimal for this type of study, such as molecular dynamics simulations, molecular docking and Monte Carlo simulations. However, this type of study is not enough for the development of biosensors.

It is of extreme importance to study the interaction of the bioreceptor with the transducing interface 116 . Examples of this are adsorption studies to support experimental methods of surface deposition 117,118 , studies on the effect of immobilization on the bioreceptor activity 119,120 or studies on the effect of the biorecognition process on the physio-chemical environment and, consequently, the transducer [START_REF] Choodet | Combined in Silico and in Vitro Study of an Aptasensor Based on Citrate-Capped AuNPs for Naked-Eye Detection of a Critical Biomarker of Oxidative Stress[END_REF]121,122 . Unfortunately, studies on the interaction biomolecule-surface (or biomolecule-nanoparticle) are very rare to this date.

The following studies discuss the effect of the addition of spacers on a bioreceptor effectiveness published in recent years, a very important factor in the development of biosensors. In 2017, Thyparambil et al 123 reported a study on the affinity of a small peptide to ochratoxin A by using biased exchange metadynamics coupled with Markov State Model to elucidate the effects of site-specific modification on the recognition process.

NF04, a designed mimotope that simulates a sensitive antibody to ochratoxin A 124 , was modified by adding six histidine residues on the N-and C-termini. They found that the Nmodified peptide had a lower affinity to the target by one order of magnitude while the Cmodified peptide had a lower affinity by five orders of magnitude. These results were in agreement with previous experimental results 125 . Figure 15 was extracted from their study 123 showcasing their results. In 2019, Rodriguez-Salazar et al 126 published a study on the in-silico design of a mixed peptide receptor sensitive to carboxyhemoglobin. This was done by performing molecular docking the target on the active site of an enzyme that binds to the target (cytochrome b561) and mimicking the spatial disposition of the amino acids that interact with the target directly while carefully adding methylene bridges between them, arriving at the conclusion that bulkier molecules tend to have a better affinity with the target due to the stabilization of the active site, to the point that polyethylene is added on both extremities of the new polymer. One issue of this work, however, is that the authors do not investigate the conformational space nor the folding process of the new proposed mixed polymer, instead assuming that the polymer is capable of assuming the required conformation for the recognition of the target. Nevertheless, this work is interesting from the point of view of the effects of adding molecular bulk to a bioreceptor's affinity to the target. In 2020, the same group published a similar study where they proposed a mixed peptide receptor sensitive to dopamine 127 .

The following studies discuss the interaction of DNA molecules with gold (both nanoparticles and surface). In 2018, Phanchai et al 121 first reported the folding of a DNA aptamer sensitive to ochratoxin A using molecular dynamics by simply running simulations of the aptamer with the target in solution and analyzing the average distance between the target, as well as the formation of H-bonds. The authors found the binding sites of the aptamer narrowed down to three nucleobases, as well as that the target would bind tightly to the aptamer after 20 ns of simulation, what was not observed when they ran simulations with other similar molecules, such as ochratoxin B and warfarin.This is interesting for this work as they were capable of narrowing down the locus of interaction to only three nucleobases, although in this thesis we attempted a full molecular docking study on our aptamer to find the accurate structure of the hairpin-ligand complex. The authors have then analyzed the adsorption of the same aptamer on citrate-capped gold nanoparticles and ran simulations with the nanoparticles bound to the aptamer in presence of ochratoxin A. They have found that the recognition of the target molecule by the aptamer renders the gold nanoparticles more exposed to the environment and therefore more prone to agglutination, which the authors propose as the mechanism of a colorimetric aptasensor based on the system (see Figure 16). In 2019, the same group published similar studies with an aptamer sensitive to 8-oxo-dG, using citrate-capped [START_REF] Choodet | Combined in Silico and in Vitro Study of an Aptasensor Based on Citrate-Capped AuNPs for Naked-Eye Detection of a Critical Biomarker of Oxidative Stress[END_REF] and cysteamine-capped 122 gold nanoparticles, this time with experimental validation in both cases. Even though the modalities of interaction aptamer-surface are effectively different from the one investigated in this thesis (adsorption instead of grafting), these studies propose that it is possible to study conformational switches of aptamers by molecular dynamics, as well as to model their interaction with gold structures. In 2020, Drew and Walsh 128 published a work on which they studied the adsorption of a DNA hairpin to a {111} gold surface in different temperatures and ionic forces. They have found that, while the ionic force didn't influence much the results, the temperature and presence of the gold surface played important roles in the disorder of the system: at 300 K, the hairpin retained an ordered state (measured by the number of stacking bases and number of h-bonds) both in solution and adsorbed to the surface, while at 350 K, the adsorbed hairpin showed a more disordered conformation than in solution. At 400 K, both the adsorbed and in solution hairpins are highly disordered, with a difference in general structures: the adsorbed hairpin attains a flat structure over the surface while the hairpin in solution maintained a globular structure. Examples of the final structures found by the authors for the adsorbed hairpin can be seen on Figure 17. This study is interesting to us as it uses the same type of gold surface (although not covered by a molecular monolayer) and its interactions to DNA hairpins, which are also of extreme importance in this thesis. More directly related to the biosensor topic of this thesis, in 2020, Lecot et al 118 published a work on the molecular modelling of a streptavidin/biotin biosensor. They studied the adsorption of streptavidin in silica surfaces that are functionalized with different types of silane monolayers and the effect of these surface treatments on the conformation of streptavidin and therefore the interaction between the streptavidin and biotin. This was done by steered molecular dynamics simulations, in which a force is applied to the biotin to remove it from the binding site of the streptavidin and the force-time curve was used to compare the difference in rupture forces between the different types of surface. Their results show that alkyl-silane treatments were the most efficient in keeping the streptavidin adsorbed while maintaining its functionality. The authors also found that the adsorption to the surface decreases the effectiveness of the binding between streptavidin and biotin, with rupture forces 2 to 4 times weaker than in solution, with the best results in maintaining adsorption and functionality found on the surface functionalized with nbutyltrimethoxysilane molecules.

Badhe et al 129 published in early 2021 a work on the in silico design of peptides sensitive to the receptor binding domain of SARS-CoV-2 and their interaction with said viral protein, as well as modifying the novel peptides (the potential bioreceptors) to functionalize graphene and carbon nanotubes in hopes of using these two materials as supporting structures in the biosensors. Figure 18 is extracted from this work showcasing their protocol for the design of the SARS-CoV-2 biosensor.

Molecular Mechanics

Computational Chemistry is the domain of chemistry on which we use When it comes to quantum mechanics, there are two different ways of calculation of the chemical properties of a system: ab-initio calculations and semiempirical molecular orbital methods. These methods solve numerically the Schrödinger equations, taking into consideration the electrons of the system.

Obviously, these methods, while are expected to be the most precise, can quickly become unviable as we increase the size of our system of interest and the time frame we want to study.

In order to decrease the computational costs while studying larger systems for longer time frames, semi-empirical methods simplify the computation in many different ways, such as only dealing with explicit valence electrons or data-fitting to experimental data. Both ab-initio and semi-empirical methods are commonly used to study phenomena taking place in the electronic scale, such as the formation and breaking of bonds or electronic rearrangements, with ab-initio methods being able to study systems with tens of atoms over time scales in picoseconds and semi-empirical methods being able to study a few hundreds of atoms in up to a nanosecond time scale. Among quantum mechanics methods, one traditional ab initio method is the Hartree-Fock method, that follows a series of approximations to simplify the numerical solving of the Schrödinger equation. Today, the most employed is the Density Functional Theory (DFT), which is widely applied for studying the electronic structure of molecules and crystal lattices. This method is based on the assumption that the energy of an electronic system can be described as a functional of its density.

When studying larger systems in longer time scales, the preferred approach for defining the interactions between particles is molecular mechanics (MM). In molecular mechanics, we assume the validity of the Born-Oppenheimer approximation, which dictates that the wave functions of the electrons and the atomic nuclei can be dealt with separately due to the enormous difference in mass between the nuclei and the electrons, so the electrons are not taken into account and charge is defined directly at the atom centers. All interatomic interactions are defined by appropriate potential energy functions that are collectively called a force field. The potentials are defined based on experimental data or localized ab-initio calculations and are applied using classical mechanics. Molecular mechanics is mostly used to study processes involving a larger number of atoms or molecules, such as surface interactions, dynamics of biomolecules and phase transitions.

General concepts about force fields

Basically, a force field consists of a series of potential energy functions that describe all the interactions between all atoms in a system, either through bonding energy (relative to the chemical bonds) or non-bonding energy. Therefore, the functional expressing the total energy in a classic force field can be described as the following:

𝐸 𝑡𝑜𝑡𝑎𝑙 = 𝐸 𝑏𝑜𝑛𝑑𝑖𝑛𝑔 + 𝐸 𝑛𝑜𝑛-𝑏𝑜𝑛𝑑𝑖𝑛𝑔 (1)
𝐸 𝑏𝑜𝑛𝑑𝑖𝑛𝑔 and 𝐸 𝑛𝑜𝑛-𝑏𝑜𝑛𝑑𝑖𝑛𝑔 can both be further developed in a sequence of individual components. 𝐸 𝑏𝑜𝑛𝑑𝑖𝑛𝑔 can be divided into the bond stretching energy (𝐸 𝑏𝑜𝑛𝑑 ), the angular bending energy (𝐸 𝑎𝑛𝑔𝑙𝑒 ), dihedral torsion energy (𝐸 𝑑𝑖ℎ𝑒𝑑𝑟𝑎𝑙 ) and improper torsion (𝐸 𝑖𝑚𝑝𝑟𝑜𝑝𝑒𝑟 ). As a side note, 𝐸 𝑑𝑖ℎ𝑒𝑑𝑟𝑎𝑙 deals with the energy involved in the torsion of a group of 4 atoms successively bonded i-j-k-l, while 𝐸 𝑖𝑚𝑝𝑟𝑜𝑝𝑒𝑟 deals with the torsion of 4 non-successively bonded atoms, the latter being useful for imposing planar structures. All these components take into account the deviation of the system from a reference structural value and usually take form in one or multiple oscillator-type potential wells. 𝐸 𝑛𝑜𝑛-𝑏𝑜𝑛𝑑𝑖𝑛𝑔 can be divided in two main components: van der Waals interactions (𝐸 𝑣𝑑𝑊 ) and electrostatic interactions (𝐸 𝑒𝑙𝑒𝑐 ). Certain force fields also add a third component for H-bond interactions ( 𝐸 ℎ𝑏𝑜𝑛𝑑 ). All these components are schematized on Figure 1. The final functional for the total energy of the system is therefore defined as:

𝐸 𝑡𝑜𝑡𝑎𝑙 = 𝐸 𝑏𝑜𝑛𝑑 + 𝐸 𝑎𝑛𝑔𝑙𝑒 + 𝐸 𝑑𝑖ℎ𝑒𝑑𝑟𝑎𝑙 + 𝐸 𝑖𝑚𝑝𝑟𝑜𝑝𝑒𝑟 + 𝐸 𝑣𝑑𝑊 + 𝐸 𝑒𝑙𝑒𝑐 + 𝐸 ℎ𝑏𝑜𝑛𝑑 (2)
The well-functioning of a force field heavily depends on the adequate parametrization of each of these components. Parameters can be obtained from abinitio calculations, experimental data and the fitting of simulation towards experimental results. Examples of typical parameters include atomic masses, van der Waals radii, average bond lengths, partial atomic charges, among others.

Another distinction between force fields is whether they are polarizable or not.

Polarizable force fields take into account the electrostatic environment of the atom when defining its effective charge at a certain moment. When compared to classical force fields, polarizable force fields are particularly interesting when we work with heavier and more polarizable atoms, however one needs to take into account that the repeated recalculation of the electrostatic field in the system may increase the computational costs of the simulation up to the point it becomes unviable.

In our work, we have used the classical Amber force fields OL15 1 for the DNA structures and the GAFF2 2 generic force field for the organic molecules (MCH and AMP). Other examples of force fields available are CHARMM 3 , OPLS 4 , MMFF 5 and UFF 6 .

The Amber force fields

The Amber force fields were first developed by Peter Kollman's group for the simulation of biomolecules 2,7 , such as proteins and nucleic acids. The validity of these force fields has been widely proven in literature [8][9][10] . The functional of the potential energy is composed of the following terms: bond stretching, angle stretching, torsions, van der Waals and electrostatic interactions. The term relative to the bond stretching is the following:

𝐸 𝑏𝑜𝑛𝑑 = ∑ 𝑘 𝑖 2 (𝑙 𝑖 -𝑙 𝑖0 ) 2 𝑏𝑜𝑛𝑑𝑠 ( 3
)
This term is an application of Hooke's law for a harmonic potential, therefore considering the chemical bonds to be equivalent to springs on which we attach one atom at each extremity. For a bond i, the parameter 𝑘 𝑖 is the force constant of the bondstretching force constant, 𝑙 𝑖 is the length of the bond at a given moment and 𝑙 𝑖0 is the reference length of the bond on which the potential energy is at the minimum.

The angle stretching term is also a harmonic potential of similar structure. This potential is given by the following equation:

𝐸 𝑎𝑛𝑔𝑙𝑒 = ∑ 𝑎 𝑖 2 (𝜃 𝑖 -𝜃 𝑖0 ) 2 𝑎𝑛𝑔𝑙𝑒𝑠 (4)
In this equation, for an angle i formed between two consecutive bonds at a given atom, the parameter 𝑎 𝑖 is the force constant associated with the harmonic potential, while 𝜃 𝑖 is the angle at a given moment and 𝜃 𝑖0 is the reference angle on which the potential energy relative to this angle is at its minimum.

For the torsion energies, the equation is the following:

𝐸 𝑑𝑖ℎ𝑒𝑑𝑟𝑎𝑙 = ∑ 𝑉 𝑛 2 (1 + 𝑐𝑜𝑠(𝑛𝜔 -𝛾)) 𝑡𝑜𝑟𝑠𝑖𝑜𝑛𝑠 (5)
In this equation, 𝑉 𝑛 is a parameter referred to as "energy barrier height" of the torsion and it gives an indication of the energy barrier present in the rotation. The parameter 𝑛 is called multiplicity, and it is the number of points on which the function attains a minimum as we rotate it 360°. The value 𝜔 represents the torsion angle at a given moment and 𝛾 is the phase factor, which is the value of a torsion angle on which the potential energy is at its minimum. Certain torsions, such as the O-C-C-O torsion have their potential defined by the combination of two terms in order to better model the differences in energy of the minima. Eventually other terms are added to take into account certain planar structures in the system.

The final term in the functional that references the non-bonding interactions, that is, both the van der Waals and the electrostatic interactions and it is as given: When it comes to the electrostatic interactions, the parameters 𝑞 𝑖 and 𝑞 𝑗 are the atomic charges while 𝜀 0 is the vacuum permittivity. As a side note, hydrogen bonds are no longer explicitly taken into account by the Amber force fields, but in other force fields or scoring functions such as the one used by Autodock 11 for the molecular docking. It is usually represented by a 12-10 LJ potential coupled with a factor taking into account the angle formed between the donor atom, the hydrogen atom and the acceptor atom.

𝐸 𝑛𝑜𝑛-𝑏𝑜𝑛𝑑𝑖𝑛𝑔 = ∑ ∑ (𝜀 𝑖𝑗 [( 𝑅 𝑖𝑗 0 𝑟 𝑖𝑗 ) 12 -2 ( 𝑅 𝑖𝑗 0 𝑟 𝑖𝑗 ) 6 ] + 𝑞 𝑖 𝑞 𝑗 4𝜋𝜀 0 𝑟 𝑖𝑗 ) 𝑁 𝑗=𝑖+1 𝑁 𝑖=1 (
When we sum up all the components of the functional, we obtain the total potential function of the force field as the following:

𝐸 𝑡𝑜𝑡𝑎𝑙 = ∑ 𝑘 𝑖 2 (𝑙 𝑖 -𝑙 𝑖0 ) 2 𝑏𝑜𝑛𝑑𝑠 + ∑ 𝑎 𝑖 2 (𝜃 𝑖 -𝜃 𝑖0 ) 2 𝑎𝑛𝑔𝑙𝑒𝑠 + ∑ 𝑉 𝑛 2 (1 + 𝑐𝑜𝑠(𝑛𝜔 -𝛾)) 𝑡𝑜𝑟𝑠𝑖𝑜𝑛𝑠 + ∑ ∑ (4𝜀 𝑖𝑗 [( 𝜎 𝑖𝑗 𝑟 𝑖𝑗 ) 12 -( 𝜎 𝑖𝑗 𝑟 𝑖𝑗 ) 6 ] + 𝑞 𝑖 𝑞 𝑗 4𝜋𝜀 0 𝑟 𝑖𝑗 ) 𝑁 𝑗=𝑖+1 𝑁 𝑖=1 (7)

Specific force-field parameters for gold and MCH

The structural details of the gold surface and the MCH monolayer used (shown on Figure 2) were obtained from previous ab-initio studies from our team 12,13 . As we did not desire any major deformations of this optimized structure, we restrained the movement of the Au, S and C atoms in this structure during our simulations, leaving only the terminal -OH group from the MCH molecule mobile, as this is the only group in direct contact with the rest of the system. The parameters for the Au atoms were manually added to the force field, and these atoms only interact with other atoms of the system by van der Waals interactions.

We have used the AM1-BCC method 14,15 to determine the partial atom charges of the MCH molecule, with the help of the antechamber software from the Amber suite. This method consists of the combination of a semi-empirical calculation using the AM1 semi-empirical quantum calculation with bond charge correction (BCC) in order to fit the charges obtained from AM1 to a training set of molecules calculated by the Hartree-Fock method using the 6-31G* basis set.

Molecular docking

Molecular docking is a method which aims to predict host/guest association 16 .

Molecular docking permits us to systematically explore the conformational space of a ligand inside a cavity and it is a very important tool for the prediction of binding conformations and orientations of small molecule drug candidates to their targets, so one can study the affinity and the activity of this small molecule. The docking methods are usually divided in two steps: a docking algorithm to explore the conformational space of the small molecule onto the target surface and a scoring function to evaluate the interaction between the small molecule and the target and predict the biological activity of this said molecule.

Therefore, molecular docking is a very important tool in computer-aided drug design. One can use molecular docking to perform very thorough screenings of immense compound libraries, evaluate the results using the scoring functions and propose conformational possibilities on the binding between the ligand and the target.

As such methods can greatly reduce the costs (both money and timewise) of drug design in the pharmaceutical industry, there is a considerable ongoing effort towards improving the methods used in molecular docking.

There are many commercially and academically available molecular docking programs, such as DOCK 17 , AutoDock 11,18 , LigandFit 19 and rDock 20 . In our work, we have used the AutoDock software to determine the structure of the aptamer hairpin-AMP complex and use this structure as basis for the construction of the complete biosensor model. AutoDock proposes a method for the binding energy evaluation by placing the receptor on a grid and calculating the energy values at each grid point based on atom types, charges and placement relative to the grid, as schematized on Figure 3. This renders the receptor's energy contributions completely described by the grid and only the grid needs to be consulted during the simulations and scoring process. One major disadvantage of this method, however, is that the receptor is now locked in place by the grid; moving the receptor demands a recalculation of all energetic contributions on the grid. The placement of the ligand on the receptor can be done by different algorithms, some of which we will discuss in the next section.

Searching algorithm

Amongst the placement algorithms available, we will discuss the Genetic Algorithm method and its improved version, the Lamarckian Genetic Algorithm.

In the Genetic Algorithm method, we consider the degrees of freedom of the ligand (orientation, translation and conformation) as genes, called state variables. The combination of the state variables is the equivalent of the genotype of the state while the atomic coordinates are the equivalent of the phenotype. The total interaction energy of the ligand with the target is the fitness of the state, calculated by the energy function. Pairs of states are randomly mated, that is, they undergo a process of crossover, on which the descendant states inherit state variables from either parent, on top of undergoing eventual random mutations. There is a selection based on the fitness of the offspring states in a generation, and the best fit ones are then retained for the next mating step. This iteration process is continued until a determined termination criteria is met. Genetic algorithm is efficient in the global search of optimal placements and it is a more adequate choice of algorithm over other more random methods such as Monte Carlo Simulated Annealing when the ligand has many degrees of freedom, as these degrees would cause a combinatorial explosion.

The genetic algorithm can be further improved by adding a local search method step in the mating process 18 . This improved method is called Lamarckian Genetic Algorithm. This makes reference to the Lamarckian theory of evolution that assumes that certain phenotypic characteristics obtained through an individual lifetime can be inherited by the offspring of said individual. In this improved algorithm, there is an inverse mapping function that provides a genotype from a given phenotype. The Lamarckian genetic algorithm works similarly to the genetic algorithm, however, now each conformation undergoes an energy minimization before the scoring calculation, and these minimized structures assigned to their genotype by the inverse mapping function and they are the ones carried forward on the process. This accelerates even further the searching step compared to the original genetic algorithm, allowing the processing of ligands with even more degrees of freedom.

Scoring function

The scoring process in Autodock follows the subsequent function of five terms:

𝛥𝐺 = 𝐶 𝑣𝑑𝑊 ∑ (
In this equation, considering all atoms in the ligand, i, and all the atoms in the receptor, j, all five 𝐶 terms on the right side are coefficients obtained empirically from linear regression analysis of a set of complexes with known coefficients 21 . All the contributions in vacuum (van der Waals, h-bonds, electrostatic and torsion) are calculated in MM terms: a 12-6 LJ summation for van der Waals interactions, another one coupled to a factor that is a function of the angle (cos 𝜃) 4 to the h-bonds, a Coulombic potential for the electrostatic interactions and the torsion contribution is directly proportional to the number of torsions (𝑁 𝑡𝑜𝑟𝑠𝑖𝑜𝑛 ) in the ligand. The solvation component is a summation of the solvent entropy changes at solute-solvent interfaces, a modelling of the de-solvation of the sites upon binding and the hydrophobic effect.

The final score is an approximation of the binding free energy; therefore, technically the most stable structures have the lowest score, although it is common to say these are the highest scoring ones in common speech. This method has proven its efficiency for protein/ligand prediction and many successful references can be found 22 . However, only few scoring functions have been specifically designed to nucleic-acids structures so that our team developed, some years ago, a specific scoring function within Autodock for ligand/nucleic-acids with the help of artificial intelligence 23 which is still used nowadays 24 .

Molecular Dynamics

Molecular dynamics is the domain on which we numerically solve the equations found by molecular mechanics allowing the study of the evolution of our system over time. When needed, the upper limits of system size and time scales can be further pushed by grouping the atoms into larger pseudo-atoms, decreasing further the degrees of freedom of the system. This is called a coarse-grain simulation and the trade-off of this approach is the loss of information on an atomic scale; however, this method can be useful for studying even larger systems, such as large protein complexes, polymers, membranes or DNA supercoiling. Due to all the approximations involved, one needs to be especially attentive to the parametrization of coarse-grained methods and this parametrization needs to be done empirically. Due to the size of our system, we will work solely using molecular mechanics coupled with molecular dynamics to study its properties.

The application of the potential energy functions retrieved by molecular mechanics on a system and the study of its evolution over time is what we call Molecular Dynamics (MD). One can say that the domain of MD simulations was first introduced by the seminal work of Alder and Wainwright published in 1957 on which they studied phase transitions of hard spheres undergoing elastic collisions 25 . In 1960, Gibson et al discussed the radiation damage events on a crystalline copper structure using a Born-Mayer potential energy model for all repulsion forces while applying a constant force at all atoms in the extremities of the crystal towards the center of the structure to simulate the cohesive forces 26 . Later in 1964, Rahman published a work on which he studied a system of 864 atoms to simulate liquid argon using a LJ potential with resulting radial distribution functions and self-diffusion constant that correlated nicely to the experimental data 27 . In 1977, the first MD study on a protein was published by McCammon et al, where they studied the dynamics of the bovine pancreatic trypsin inhibitor in vacuum for 8.8 ps 28 . Nowadays, MD simulations is a well-established method of computational studies and in this section, we will discuss further some technical aspects regarding MD simulations.

Equations and time steps

In MD simulations, the movement of all particles is governed by Newtonian mechanics. The position of each atom i is defined by a vector 𝑟 𝑖 ⃗ ⃗ . All numerical

calculations
𝑎 𝑖 ⃗⃗⃗ = 𝐹 𝑖 𝑚 𝑖 ⃗⃗⃗ = - ∇ ⃗ ⃗ 𝐸 𝑡𝑜𝑡𝑎𝑙 𝑚 𝑖 (10)
In order to continuously analyze the system's evolution over time, we also need to keep track of each atom's position and velocity over time. In this work, the algorithm used is the leapfrog integration method 29 . In this method, both position and velocity are updated in alternating points in time every half time step. The equations used for such are the following:

𝑟 𝑖 ⃗ ⃗ (𝑡 + 𝛥𝑡) = 𝑟 𝑖 ⃗ ⃗ (𝑡) + 𝑣 𝑖 ⃗⃗ ⃗ (𝑡 + 𝛥𝑡 2 ) 𝛥𝑡 (11)
𝑣 𝑖 ⃗⃗ ⃗ (𝑡 + 𝛥𝑡 2 ) = 𝑣 𝑖 ⃗⃗⃗ (𝑡 - 𝛥𝑡 2 ) + 𝑎 𝑖 ⃗⃗⃗ (𝑡) 𝛥𝑡 (12)
where 𝑣 𝑖 ⃗⃗⃗ and 𝑟 𝑖 ⃗ ⃗ are the velocity and position of an atom i respectively, and 𝛥𝑡 is the time step of the simulation. Initial velocities are assigned to the atoms following a

Maxwell-Boltzmann distribution. The leapfrog method is also advantageous by being time-reversible and mostly conserving the total energy of the system.

When it comes to the time step for a MD simulation, 1 fs is the safest value of choice 30 , due to the usual timeframe in which relevant interactions happen in this scale.

However, there are certain algorithms such as the SHAKE 31 method that facilitate the increase of the time step to 2 fs. The SHAKE method takes into account that the bonds with higher stretch frequency are the ones involving hydrogen atoms. A constraint is added to these bonds and the position of the hydrogen atoms is therefore dependent on the position of the atoms to which they are bonded. This method is largely used in MD simulations of biomolecules due to the sheer size of the systems. Usually, the information lost during its application, that is, the vibration of the hydrogen atoms, is not relevant to these simulations, but it is important to keep in mind that information such as the infrared vibrational frequencies of these bonds cannot be obtained, and QM/MM methods are not advised for that purpose.

Free energy determination: the Umbrella Sampling technique

Getting the value of the free energy is widely considered as the most important quantity that can be attained through molecular simulation. Indeed, the free energy difference represents the driving force of any process, such as any chemical reaction, and encompasses entropy which are somewhat seen as elusive.

At first sight, as molecular dynamic computes both energy and Newton's equations, we could expect to generate a free energy curve along one reaction coordinate solely from one molecular dynamic trajectory. This requires the analyses of probability of states and thus their respective populations. However, the barrier energy between states is many times larger than the kBT energy that is accessible during a molecular dynamic simulation which, in fact, samples the phase space around the local minima and remain stagnated into a valley of potential. In reality, molecular dynamics simulations should be undertaken for trajectories on second to hour timescale instead of hundred nanoseconds, which is, to date, the classical timescale. This is obviously not tractable and so a smarter approach must be used.

According to the challenge of getting free energies, many computational techniques have been and are still developed. All of them try to enhance the sampling along the reaction coordinate. One can cite methods which:

o Prevent the exploration of "already seen states" by adding a restraint. This is the case of metadynamic 32 techniques.

o Guide the exploration along a reaction coordinate by adding a restraint which evolves along the time. The steered molecular dynamics 33 (SMD) and adaptive steered molecular dynamics (aSMD) fall in this category.

o Explore, methodically, each state defined along the reaction coordinate by making numerous molecular dynamics simulations (called windows) trapped around a value of the coordinate. This is called the umbrella sampling 34 (US) method.

In 2018, Wei and coworkers review 35 these three above-mentioned methods for interfacial systems such as biomolecules onto inorganic surfaces. Even though umbrella sampling is the one which requires more computational investigations, this technique appears to be the most reliable one 35,36 . We thus decided to use this one in this work.

One key feature of umbrella sampling is to make numerous molecular dynamics simulations, called windows, with a fixed restraints along the reaction coordinate. A harmonic restraint is then applied, giving the bias potential 𝑊(𝑡) according to the following equation:

𝑊(𝑡) = 𝑘 𝑤 (𝜉(𝑡) -𝜉 0 ) 2 (13)
where 𝜉 0 is the reference reaction coordinate whereas 𝜉(𝑡) is its value at time 𝑡. 𝑘 𝑤 is a force constant and its value could be optimized when such simulations are realized.

Umbrella sampling follows the reaction coordinate and explore this one through a series of independent MD simulations (windows). For each of these windows, a harmonic restraint is applied to force the reaction coordinate to remain close to the center of each window. From this, a series of histograms showing the sampling of the reaction coordinate values can be created. If these windows overlap, the weighted histogram analysis method (WHAM 37 ) can be employed to post process the results to remove the biasing of the restraints and recover the potential mean force (PMF), that is, the image of the free energy surface following the reaction coordinate.

In a way, the US technique is merely a juxtaposition of many independent molecular dynamic simulations with a harmonic restraint. For each simulation, the value for the minimum of the harmonic restraint potential evolves along the reaction coordinate to consider the whole process. The heating phase is made in the NTP ensemble whereas the equilibration and production are in the NVT ensemble. A schematic of this protocol is shown on Figure 4.

Water models for explicit solvation

One approach to model a solvent is to represent it explicitly in an all-atom model inside the system to be studied. However, as a drastic number of solvent molecules are added, it is of upmost importance to use the most appropriate potential model to describe the intermolecular interactions involving the solvent 38 . There are many models for the water molecule that try to solve the question of the potential energy, such as SPC 39 , SPC/E 40 , TIP3P 41 and TIP4P 41 .

One main feature of all these water molecules is the usage of a LJ sphere centered at the oxygen atom and positive charges situated at both hydrogen atoms.

Models on which the negative charge is located directly at the oxygen are called 3-site models, it is the case of the SPC, SPC/E and TIP3P models for example. Certain models choose to place their negative charges in a dummy atom at the bisector of the angle HOH, these called 4-site models and TIP4P is an example of such model. Some other models go even further and add a dummy atom at each of the lone pairs of the oxygen to simulate their charges (5-site models) and others combine both the bisector approach and the lone pairs approach (6-site models). These are all illustrated in Figure 5. In the simulations presented in this work, the model of choice was TIP3P. This model is very often used in the molecular dynamics simulations of biomolecules, with good results. It is also a good compromise between computational costs and energy results, as models with more sites imply more interactions to be computed per molecule, which when multiplied by the number of solvent molecules in the system can increase computational times significantly. The parameters 41 of the TIP3P model are as following:

Periodic boundary conditions

When running MD simulations, one needs to take into account the effects of boundary on the atoms. If an atom is at the boundary of a finite simulation space, it will only undergo forces that originate from the simulation space when compared to an atom in the bulk of the system, undergoing forces from all directions (Figure 6).

Therefore, the atoms are not undergoing the same forces and the system's isotropy is broken. One way to solve this issue is to extend this system indefinitely. However, it is impracticable to simulate an infinite system without any simplification. For simulations on which we explicitly add the solvent in the atomic level, this simplification is achieved by what we call Periodic Boundary Conditions (PBC). In PBC, the simulated system is replicated in all directions forming an infinite lattice. Anytime an atom leaves the initial system towards a neighboring replica, an equivalent atom is added to the initial system entering it from the opposite boundary, as represented on Figure 7. As all replicas undergo the same processes simultaneously, this ensures that all simulation boxes always have the same number of atoms at all times. The PBC method isn't without any complications since the precise calculation of all the forces acting on an atom must consider the atoms of all neighboring replicas, which now are infinite. To overcome this problem, we can add a cutoff distance for all explicit interactions to truncate the potential energy calculation. That is, to an arbitrary atom in the system, we only explicitly calculate the interactions for the atoms that are closer than a defined cutoff distance. This is done because the atoms closer to this arbitrary atom are the ones that impact the most the potential energy and the force applied on the arbitrary atom at any time. The interactions with the atoms further from the cutoff distance are simplified as a homogeneous medium. Long range corrections of the electrostatic interactions of this method are introduced by the Particle Mesh Ewald method 42 . One limitation of the Particle Mesh Ewald that is relevant to this work, however, is that the system needs to be electrically neutral, so that counter ions should be added in order for ionic systems to reach neutrality.

Implicit solvation

Implicit solvation aims to model the solvent as a continuous medium instead of explicitly represent water molecules. One important factor to take into account when working with implicit solvation is the solvation free energy, that is, the energy required to move a molecule from the gas phase (effectively surrounded by vacuum) to a solvated phase. This energy is important, for example to estimate the binding free energy involved in the formation of a complex, as exemplified on the thermodynamic cycle shown on Figure 8. The solvation free energy ( 𝛥𝐺 𝑠𝑜𝑙 ) can be decomposed in three main components: 𝛥𝐺 𝑒𝑙𝑒𝑐 , for the electrostatic interactions between solvent and solute, 𝛥𝐺 𝑣𝑑𝑊 for the van der Waals interactions between solvent and solute and 𝛥𝐺 𝑐𝑎𝑣 , which is the energy required to form a cavity to contain the solute inside the solvent. The latter also includes an entropic factor connected to the reorganization of the solvent molecules to accommodate the solute as well as the pression from solvent during the cavity creation. If the solvent and the solute form a hydrogen bond between each other, a fourth term 𝛥𝐺 ℎ𝑏𝑜𝑛𝑑 may be added to 𝛥𝐺 𝑠𝑜𝑙 .

𝛥𝐺 𝑠𝑜𝑙 = 𝛥𝐺 𝑒𝑙𝑒𝑐 + 𝛥𝐺 𝑣𝑑𝑊 + 𝛥𝐺 𝑐𝑎𝑣 (14)

There are two main methods to define 𝛥𝐺 𝑒𝑙𝑒𝑐 : Generalized Born (GB) models and the Poisson-Boltzmann (PB) models. These methods are detailed below.

Generalized Born (GB) model

The GB model is considered accurate and has a low computational cost 43 . This model considers all the atoms in a molecule with 𝑁 atoms as a sphere of radius 𝐴 𝑖 and charge 𝑞 𝑖 . Still considering the solvent as a homogeneous medium of relative permittivity 𝜀, the model describes the total electrostatic free energy 𝐺 𝑒𝑙𝑒𝑐 as:

𝐺 𝑒𝑙𝑒𝑐 = ∑ ∑ 𝑞 𝑖 𝑞 𝑗 𝜀𝑟 𝑖𝑗 - 1 2 (1 - 1 𝜀 ) 𝑁 𝑗=𝑖+1 𝑁 𝑖=1 ∑ 𝑞 𝑖 2 𝐴 𝑖 𝑁 𝑖=1 (15)
𝑟 𝑖𝑗 being the distance between atoms i and j. However, the first term of this equation can be further developed as follows:

∑ ∑ 𝑞 𝑖 𝑞 𝑗 𝜀𝑟 𝑖𝑗 𝑁 𝑗=𝑖+1 𝑁 𝑖=1 = ∑ ∑ 𝑞 𝑖 𝑞 𝑗 𝑟 𝑖𝑗 𝑁 𝑗=𝑖+1 𝑁 𝑖=1 -(1 - 1 𝜀 ) ∑ ∑ 𝑞 𝑖 𝑞 𝑗 𝑟 𝑖𝑗 𝑁 𝑗=𝑖+1 𝑁 𝑖=1 (16)
That when we substitute back to the 𝐺 𝑒𝑙𝑒𝑐 equation, we obtain:

𝐺 𝑒𝑙𝑒𝑐 = ∑ ∑ 𝑞 𝑖 𝑞 𝑗 𝑟 𝑖𝑗 𝑁 𝑗=𝑖+1 𝑁 𝑖=1 -(1 - 1 𝜀 ) ∑ ∑ 𝑞 𝑖 𝑞 𝑗 𝑟 𝑖𝑗 𝑁 𝑗=𝑖+1 𝑁 𝑖=1 - 1 2 (1 - 1 𝜀 ) ∑ 𝑞 𝑖 2 𝐴 𝑖 𝑁 𝑖=1 (17)
One important detail to notice, however, is that the first term of this equation is equal to the Coulomb interaction of two charges in vacuum, so we can obtain the following:

𝛥𝐺 𝑒𝑙𝑒𝑐 = 𝐺 𝑒𝑙𝑒𝑐 -∑ ∑ 𝑞 𝑖 𝑞 𝑗 𝑟 𝑖𝑗 𝑁 𝑗=𝑖+1 𝑁 𝑖=1 = -(1 - 1 𝜀 ) ∑ ∑ 𝑞 𝑖 𝑞 𝑗 𝑟 𝑖𝑗 𝑁 𝑗=𝑖+1 𝑁 𝑖=1 - 1 2 (1 - 1 𝜀 ) ∑ 𝑞 𝑖 2 𝐴 𝑖 𝑁 𝑖=1 (18)
When using the Amber software, however, this equation is approximated to the following:

𝛥𝐺 𝑒𝑙𝑒𝑐 = - 1 2 ∑ ∑ 𝑞 𝑖 𝑞 𝑗 𝑓 𝐺𝐵 (𝑟 𝑖𝑗 ,𝑅 𝑖 ,𝑅 𝑗 ) (1 - 𝑒 -𝜅𝑓 𝐺𝐵 𝜀 𝑁 𝑗=𝑖+1 𝑁 𝑖=1 ) (19)
On which 𝜅 is the Debye-Hückel screening parameter ( 𝜅 = 0 for solvents without ions), 𝑅 𝑖 and 𝑅 𝑗 are the effective Born radii of atoms i and j, and 𝑓 𝐺𝐵 is a smoothing function to render the 𝛥𝐺 𝑒𝑙𝑒𝑐 well-behaved and differentiable in respect to 𝑟 𝑖𝑗 from close range when 𝑟 𝑖𝑗 tends to 0, to when the atoms are far enough from each other to be considered only as point charges by Coulomb's law. The function 𝑓 𝐺𝐵 is defined as the following:

𝑓 𝐺𝐵 = √ 𝑟 𝑖𝑗 2 + 𝑅 𝑖 𝑅 𝑗 𝑒 -𝑟 𝑖𝑗 2 4𝑅 𝑖 𝑅 𝑗 (20)
Finally, we need to define the effective Born radii of the atoms. The effective Born radius of an atom can be thought of as the distance between the atom and the surface of the solvent cavity. It is defined as the radius that when inserted to the formulae above, it would give the same 𝛥𝐺 𝑒𝑙𝑒𝑐 value as a hypothetical molecule having the same dielectric boundaries but having only a single charge 𝑞 𝑖 at the same position as the original 44 . This radius is usually estimated by solving an integral of the Coulomb field over the volume of the molecule 44 .

Poisson-Boltzmann (PB) method

The Poisson-Boltzmann (PB) method uses the Poisson equation to relate the electrostatic potential 𝜙(𝑟) in a medium of permittivity 𝜀 and charge density 𝜌(𝑟) as the following 45 :

∇ 2 𝜙(𝑟) = -4𝜋𝜌(𝑟) 𝜀 (21)
In the case on which we have mobile ions of charge 𝑧 𝑖 in the solution, we have to take into account the movement of these ions under the influence of 𝜙(𝑟), as this movement influences 𝜙(𝑟) itself. We assume that these ions follow a Boltzmann distribution and combine both equations to form the Poisson-Boltzmann equation. This equation is as follows:

∇𝜀∇𝜙(𝑟) = -4𝜋𝜌(𝑟) -4𝜋 ∑ 𝑧 𝑖 𝑐 𝑖 𝑒 -𝑧 𝑖 𝜙(𝑟) 𝑘 𝐵 𝑇 𝑖 ( 22
)
being 𝑐 𝑖 the concentration of ion i, 𝑘 𝐵 the Boltzmann constant and 𝑇 the temperature.

In this model, we consider the solute as a group of point charges and we calculate the interactions between the solute and the solvent using the PB electrostatic theory, which unfortunately has no analytical solution. Therefore, to numerically solve it, one can use for example a finite-difference algorithm in a cubic grid 46 . The APBS software 47 , for example, places a cubic grid over the solute and its surroundings and assigns values of electrostatic potential, charge density, permittivity and ionic strength on each point of the grid.

One essential step of this calculation is the allocation of proper permittivity values to each point in the grid; therefore, it is of utmost importance to properly attribute which points are inside the solute and which are in the solvent. The definition of the solute-solvent boundary can be done in two ways: molecular surface or accessible surface. To the grid points outside this surface, it is allotted the permittivity of the solvent 𝜀 𝑒𝑥 (around 80 for water), while to grid points inside the surface it is allotted the permittivity of the solute molecule 𝜀 𝑖𝑛 (usually between 2 and 4).

Once we have this grid settled, we can use the PB method to calculate the 𝛥𝐺 𝑒𝑙𝑒𝑐 relative to many different processes such as solvation and the formation of intermolecular complexes. When it comes to calculating the 𝛥𝐺 𝑒𝑙𝑒𝑐 in a solvation process, we can consider that the electrostatic energy of a particle of charge 𝑞 𝑖 in a potential 𝜙 𝑖 is equal to 𝑞 𝑖 𝜙 𝑖 . Therefore, 𝛥𝐺 𝑒𝑙𝑒𝑐 can be obtained by the following:

𝛥𝐺 𝑒𝑙𝑒𝑐 = 1 2 ∑ 𝑞 𝑖 𝑖 (𝜙 𝑖 𝑠𝑜𝑙𝑣𝑒𝑛𝑡 -𝜙 𝑖 𝑣𝑎𝑐𝑢𝑢𝑚 ) (23)
The Amber software approximates this equation using what is called the Analytical Linearized Poisson-Boltzmann equation (ALPB) 48,49 , which is the following:

𝛥𝐺 𝑒𝑙𝑒𝑐 ≈ 𝛥𝐺 𝐴𝐿𝑃𝐵 = - 1 2 ( 1 𝜀 𝑖𝑛 - 1 𝜀 𝑒𝑥 ) 1 1+𝛼𝛽 ∑ 𝑞 𝑖 𝑞 𝑗 ( 𝑖𝑗 1 𝑓 𝐺𝐵 - 𝛼𝛽 𝐴 ) (24)
on which 𝛼 is a constant of value 0.571412, 𝛽 is the value of the ratio 𝜀 𝑖𝑛 𝜀 𝑒𝑥 ⁄ and 𝐴 is the effective electrostatic size of the molecule. This value is related to the effective Born radius by the following relation:

𝐴 = lim 𝜀 𝑖𝑛 →∞ 𝑅 𝑖 , 𝜀 𝑒𝑥 = 1 (25)
This method can also be used to define the electrostatic contribution to the binding free energy of a complex. This will be further discussed on section 5.2.

Non-electrostatic contributions

As a reminder, the solvation free energy 𝛥𝐺 𝑠𝑜𝑙 can be mainly decomposed in three components: 𝛥𝐺 𝑒𝑙𝑒𝑐 being the electrostatic contribution, 𝛥𝐺 𝑣𝑑𝑊 being the van der Waals contribution and 𝛥𝐺 𝑐𝑎𝑣 being the cavity formation contribution:

𝛥𝐺 𝑠𝑜𝑙 = 𝛥𝐺 𝑒𝑙𝑒𝑐 + 𝛥𝐺 𝑣𝑑𝑊 + 𝛥𝐺 𝑐𝑎𝑣 ( 26
)
𝛥𝐺 𝑐𝑎𝑣 is the contribution of free energy induced by the solvent perturbation generated by the formation of a cavity. It is always positive and includes the entropic penalty associated with the reorganization of water molecules around the solute.

𝛥𝐺 𝑐𝑎𝑣 also takes into account the work necessary against the pressure of the solvent in the cavity. 𝛥𝐺 𝑣𝑑𝑊 is the van der Waals contribution of solvent (water in our work) made around the solute. 𝛥𝐺 𝑐𝑎𝑣 and 𝛥𝐺 𝑣𝑑𝑊 represent the non-electrostatic contribution (𝛥𝐺 𝑛𝑜𝑛-𝑒𝑙𝑒𝑐 ) in implicit solvation and are gathered together. 𝛥𝐺 𝑛𝑜𝑛-𝑒𝑙𝑒𝑐 is proportional to the solvent accessible surface area of the solute (𝑆𝐴𝑆𝐴), which is obtained by a rolling ball algorithm 50 . In this algorithm, a spherical probe of a certain diameter (usually 1.4 Å as it is the approximate radius of a water molecule) is rolled over the van der Waals surface of the molecule, creating two molecular surfaces: the contact surface, which is the surface formed by the areas touched by the probe and the solvent accessible surface, which is the surface formed by the center of the probe during the process. This is illustrated by Figure 9. The equation for defining 𝛥𝐺 𝑛𝑜𝑛-𝑒𝑙𝑒𝑐 is as follows:

𝛥𝐺 𝑛𝑜𝑛-𝑒𝑙𝑒𝑐 = 𝛾 𝑆𝐴𝑆𝐴 + 𝛽 (27)
being both 𝛾 and 𝛽 constants in Amber, the Molsurf 51 module calculates 𝑆𝐴𝑆𝐴 and 𝛥𝐺 𝑐𝑎𝑣 with the values for 𝛾 = 0.00542 kcal/mol.Å 2 and 𝛽 = 0.92 kcal/mol.

MD trajectory analyses 5.1 Classical trajectory analyses

We have used the VMD software 52 to visualize our trajectories and create our figures while we have performed most of our structural analyses with the cpptraj 53 module from AmberTools 7 .

RMSD measurements were made to analyze the evolution of the coordinates over time and how far these coordinates move away from a reference set of coordinates (usually in our case, the first frame of the simulation). The mass-weighted RMSD equation is:

𝑅𝑀𝑆𝐷 = √ ∑ [𝑚 𝑖 .(𝑋 𝑖 ⃗⃗⃗⃗ -𝑋0 𝑖 ⃗⃗⃗⃗⃗⃗⃗ ) 2 𝑁 𝑖=0 ] 𝑀 (28)
in which 𝑁 is the number of atoms i involved in the measurement, 𝑚 𝑖 is the atomic mass of each atom i, 𝑋 𝑖 ⃗⃗⃗ is the coordinate vector for atom i, 𝑋0 𝑖 ⃗⃗⃗⃗⃗⃗ is the coordinate vector for atom i in the reference set of coordinates and 𝑀 is the total mass of all atoms involved in the measurement.

RMSF measurements analyze how fluctuates each fragment of our structure, by measuring the average deviation over time of each particle coordinate vector to a reference set of coordinates, in our case, the average coordinates of the analyzed simulations. The RMSF value of an atom i is obtained with the help of the following equation:

𝑅𝑀𝑆𝐹 𝑖 = √ ∑ (𝑋 𝑖 ⃗⃗⃗⃗ (𝑡 𝑗 )-𝑋0 𝑖 ⃗⃗⃗⃗⃗⃗⃗ ) 2 𝑇 𝑡 𝑗 =0 𝑇 ( 29
)
in which 𝑇 is the total time analyzed, 𝑋 𝑖 (𝑡 𝑗 ) is the coordinate vector for atom i at time 𝑡 𝑗 .

The analyses of radius of gyration (𝑅 𝑔 ) were performed to estimate how folded in themselves certain DNA structures were during simulations. They were calculated using the following equation:

𝑅 𝑔 = √ ∑ 𝑋𝑀 𝑖 ⃗⃗⃗⃗⃗⃗⃗⃗ 2 𝑚 𝑖 𝑖 ∑ 𝑚 𝑖 𝑖 ( 30
)
where 𝑋𝑀 𝑖 ⃗⃗⃗⃗⃗⃗⃗ is the coordinate vector of atom i in respect to the center of mass of the total group of analyzed atoms.

We have also performed a rough estimate of a free energy cartography using two arbitrary structural coordinates in one of our systems using a script previously written by our group 13 and the following equation:

∆𝐺(𝛷 1 , 𝛷 2 ) = -𝑘 𝐵 𝑇 . ln [ 𝑃(𝜌 𝛷 1 ,𝛷 2 ) 𝑃 𝑚𝑎𝑥 (𝜌 𝛷 1 ,𝛷 2 ) ] (31)
In this equation, 𝑘 𝐵 and 𝑇 are the Boltzmann constant and temperature (300 K) respectively. 𝑃(𝜌 𝛷 1 ,𝛷 2 ) is the probability density function of two arbitrary structural coordinates 𝛷 1 and 𝛷 2 , while 𝑃 𝑚𝑎𝑥 (𝜌 𝛷 1 ,𝛷 2 ) is the most probable state.

Binding free energy calculations of complex formation: MM-PBSA and MM-GBSA

In 1998, Kollman et al 54 published the method known as Molecular Mechanics -Poisson-Boltzmann surface area (MM-PBSA), on which they move around a thermodynamic cycle to estimate the binding free energy of a complex 𝛥𝐺 𝑏𝑖𝑛𝑑 . The value of 𝛥𝐺 𝑏𝑖𝑛𝑑 is assumed to be the sum of the gas-phase energies, solvation-free energies and the entropic factors averaged over time in a MD trajectory (see Fig 8).

The gas-phase terms are obtained through MM calculations. In order to calculate the 𝛥𝐺 𝑒𝑙𝑒𝑐 value of the formation of the complex, we calculate 𝛥𝐺 𝑒𝑙𝑒𝑐 using the PB equation for the transfer of each isolated species from the solvent to a medium of permittivity 𝜀 = 𝜀 𝑖𝑛 . The free energy involved in bringing the species together is calculated using Coulomb's law in a medium of permittivity 𝜀 = 𝜀 𝑖𝑛 and then the free energy involved in the transfer of the complex from the medium of permittivity 𝜀 = 𝜀 𝑖𝑛 to the original solvent is calculated. One intrinsic problem of this method, however, is that it assumes that both molecules have the same internal permittivity 𝜀 𝑖𝑛 .

A modification has been introduced on which the solvation energy is defined by the GB calculations 54 . The model called Molecular Mechanics -Generalized Born surface area (MM-GBSA) uses the GB approach for 𝛥𝐺 𝑒𝑙𝑒𝑐 while the non-electrostatic contributions are calculated directly from the atomic coordinates. Although one can also decompose 𝛥𝐺 𝑒𝑙𝑒𝑐 to its atomic components using MM-PBSA, this requires more costly calculation times 55 . Using MM-GBSA for this type of calculation can achieve this type of information much faster than other methods.

MM-GBSA and MM-PBSA represent a very good balance between computational efficiency and prediction accuracy 56 . MM-PBSA and MM-GBSA differ only by their electrostatic contribution of solvent. Theoretically, the PB methods are more accurate than the GB ones. However, the Poisson-Boltzmann solvers do not converge easily so that GB methods represents generally a more reliable solution 56 .

Moreover, in drug-design strategies, MM-GBSA has been found to provide better ranking capabilities 57 ) and this preference is also enhanced when studying RNA complexes 58 ). For all these reasons, we chose the MMGBSA technique as the method for the calculation of the binding energy of the aptamer-AMP complex. These calculations were performed using MD simulations of the complex in implicit solvent under the assumption that the conformations obtained for both the aptamer and AMP during the simulation are comparable in free energy to trajectories of either separated.

The calculation performed can be described by the following series of equations:

< 𝛥𝐺 𝑏𝑖𝑛𝑑 > = < 𝐺 𝑐𝑜𝑚𝑝𝑙𝑒𝑥 > -< 𝐺 𝑎𝑝𝑡𝑎𝑚𝑒𝑟 > -< 𝐺 𝐴𝑀𝑃 > (32)
< 𝐺 > = < 𝐸 𝑔𝑎𝑠 > + < 𝐺 𝑠𝑜𝑙 > -𝑇 < 𝑆 𝑔𝑎𝑠 > (33)
< 𝐸 𝑔𝑎𝑠 > = < 𝐸 𝑏𝑜𝑛𝑑 > +< 𝐸 𝑎𝑛𝑔𝑙𝑒 > +< 𝐸 𝑡𝑜𝑟𝑠𝑖𝑜𝑛 > + < 𝐸 𝑣𝑑𝑊 > +< 𝐸 𝑒𝑙𝑒𝑐 > (34)

< 𝐺 𝑠𝑜𝑙 > = < 𝐺 𝑒𝑙𝑒𝑐 > + < 𝐺 𝑛𝑜𝑛-𝑝𝑜𝑙𝑎𝑟 > (35)
-𝑇 < 𝑆 𝑔𝑎𝑠 > = -𝑇 < 𝑆 𝑣𝑖𝑏 > -𝑇 < 𝑆 𝑡𝑟𝑎𝑛𝑠 > -𝑇 < 𝑆 𝑟𝑜𝑡 > (36)
on which 𝐺 𝑐𝑜𝑚𝑝𝑙𝑒𝑥 , 𝐺 𝑎𝑝𝑡𝑎𝑚𝑒𝑟 and 𝐺 𝐴𝑀𝑃 are the free energies of the complex, the aptamer and the AMP molecule, respectively. <A> means an average of A values during the molecular dynamics trajectory. We can obtain each by the summation of the internal energy in gas-phase (𝐸 𝑔𝑎𝑠 ), the solvation free energy (𝐺 𝑠𝑜𝑙 ) and the entropic factor (𝑇𝑆 𝑔𝑎𝑠 ). The factor 𝐸 𝑔𝑎𝑠 factor is deduced from the force field potential energy functions. As we are using the "same trajectory method" (STM), the internal energy values (𝐸 𝑏𝑜𝑛𝑑 + 𝐸 𝑎𝑛𝑔𝑙𝑒 + 𝐸 𝑡𝑜𝑟𝑠𝑖𝑜𝑛 ) are identical in both the complex and the isolated aptamer and AMP molecule, so their contribution to 𝛥𝐺 𝑏𝑖𝑛𝑑 is none. The solvation free energy (𝐺 𝑠𝑜𝑙 ) is calculated as the sum of the electrostatic contribution (𝐺 𝑒𝑙𝑒𝑐 ) and the non-polar contribution (𝐺 𝑛𝑜𝑛-𝑝𝑜𝑙𝑎𝑟 ). The non-polar contribution, as discussed in the previous section, is calculated from the solvent accessible surface area (𝑆𝐴𝑆𝐴). The electric contribution is calculated by the GB method. The entropic term is decomposed in three terms: 𝑇𝑆 𝑣𝑖𝑏 for the vibrational component, 𝑇𝑆 𝑡𝑟𝑎𝑛𝑠 for the translational component and 𝑇𝑆 𝑟𝑜𝑡 for the rotational component. The vibrational component (𝑇𝑆 𝑣𝑖𝑏 ) is calculated using normal mode analysis 59 , while the latter two are direct functions of the mass and the moment of inertia, respectively.

This entire process was done with the help of the built-in mm_pbsa.pl script in Amber, which renders the whole process automatic using only one input file for the extraction and definition of the structural components from a MD trajectory in explicit solvent and another for the parametrization of the method.

Specificities of the Au/MCH system

Simulations at interfaces are not as usual as for isotropic systems, and the conditions for running them are not as straightforward. In fact, molecular mechanics is well suited for molecules but not intended for metallic surfaces. However, surface functionalization plays an enormous part in biosensors, therefore simulations at interfaces cannot be avoided. In the next sections, the methodological approaches for dealing with the interface will be discussed.

Preparation of the surface

The gold surface model used is a Au{111} surface with two atomic layers functionalized with a 6-mercapto-1-hexanol (MCH) self-assembled monolayer with a saturation coverage c(2×2). The structural formula of this molecule is shown on Figure 10. This surface was modelled based on an ab initio study previously made by the group 12,13 , which provided us a POSCAR file, that was converted to a PDB file with the help of the VMD software 52 . One important detail of this structure is that the MCH molecule does not stay in a 90° angle to the surface, due to the interaction of the free doublets of the sulfur atom with the gold surface. The visualization of this PDB file is shown below on Figure 11: We have propagated this cell using the PropPDB software from AmberTools 7 , that gives us a surface in the shape of a bent parallelogram. This surface is then manually edited to generate a PDB file containing a rectangular surface containing 2220 gold atoms and 350 MCH molecules in total, of size 9.1 x 8.9 nm 2 and grafting density of around 4.32 molecules per mm 2 , removing unwanted fragments and completing molecules that might have been cut in half in this process. This is the surface we used to build our structures in this work. This surface is shown on Figure 12. We have inserted extra force field parameters to add a strong LJ potential to keep the Au atoms in place. We have built the MCH molecules using Maestro [START_REF] Kergoat | A Water-Gate Organic Field-Effect Transistor[END_REF] to generate two mol2 files, one containing one molecule of MCH as is on the surface; that is, with a terminal -S -group and another containing the same molecule with a terminal -SH group. These two mol2 files are used as input files for the antechamber program of the Amber suit 7 to generate two files containing the partial atom charges of both species. We have fused both, using the calculated charges of the protonated species and applying to the deprotonated structure. The charges of the -SH group were combined on the sulfur atom of the deprotonated structure, as this atom is more polarizable and will be deep in the SAM structure, therefore less likely to interact with the environment. Furthermore, this atom is constrained in our MD simulations. This final file containing the partial charges of the MCH molecule was then used as input file to the parmchk program of AmberTools 7 to generate the parameters that might be missing in the GAFF2 force field 61 .

When it comes to grafting the DNA to the surface, the DNA was first prepared to be grafted to the surface by manually and slightly moving the last 3 to 4 bases of the sequence on the terminus to be grafted to render it sterically available using Maestro [START_REF] Kergoat | A Water-Gate Organic Field-Effect Transistor[END_REF] .

The grafting is done by removing the hydrogen atom of the hydroxide group of the selected MCH molecule and the hydroxide terminal group of the DNA, then manually drawing the phosphate group binding the oxygen atoms of both groups with the help of the xleap program from AmberTools 7 . The force-field parameters for the junction were adapted from the OL15 force field 1 and the partial charges were manually added and adjusted to render the entire system charge as an integer value.

Once the system was ready for solvation, we added TIP3P water 62 solvent cuboids anisotropically avoiding the addition of water molecules outside of the section of the xOy plane that contains the gold surface, as well as Na + ions up to the concentration of 160 mM and Cl -ions to counter the excessive charges, bringing the system back to a total neutral charge.

Charges on surface: the electrode simulations

We designed our biorecognition structure having in mind its possible utilization in an EGOFET-type biosensor. In such device, few hundreds of mV are applied between gate and source electrodes to switch-on the drain current. Therefore, it is important to analyze the effects of electric charges on the Au surface on the DNA conformations and by consequence, on the ionic distribution.

To find the electrical charge at the surface, we will use the hypothetical closed circuit containing the capacitor described in Figure 13 to better visualize certain aspects we need to take into account: If we consider the entire capacitor formed between the gold and the organic semiconductor (OSC) interfaces, the charges per area σ [C.cm -2 ] at the gate and OSC interfaces are equal in absolute values. This absolute value |σ| can be defined by the total capacitance of the system 𝐶 𝑡𝑜𝑡 times the absolute value of V. Being 𝐶 𝑡𝑜𝑡 the total capacitance of the system, is can also be defined by the following equation:

𝐶 𝑡𝑜𝑡 -1 = 𝐶 𝐺𝑎𝑡𝑒/𝐸𝑙𝑒𝑐𝑡𝑟𝑜𝑙𝑦𝑡𝑒 -1 + 𝐶 𝑂𝑆𝐶/𝐸𝑙𝑒𝑐𝑡𝑟𝑜𝑙𝑦𝑡𝑒 -1 (37)
being 𝐶 𝐺𝑎𝑡𝑒/𝐸𝑙𝑒𝑐𝑡𝑟𝑜𝑙𝑦𝑡𝑒 and 𝐶 𝑂𝑆𝐶/𝐸𝑙𝑒𝑐𝑡𝑟𝑜𝑙𝑦𝑡𝑒 the capacitances of the gate/electrolyte and OSC/electrolyte interfaces, respectively.

Practically, EGOFETs are usually manufactured taking into account on which of the interfaces the phenomena of interest will happen: for systems on which we are interested on effects occurring on the gate/electrolyte interface, we design the EGOFET so that the channel length and width on the OSC/electrolyte interface are much bigger, increasing 𝐶 𝑂𝑆𝐶/𝐸𝑙𝑒𝑐𝑡𝑟𝑜𝑙𝑦𝑡𝑒 to the point on which we can directly assume that only 𝐶 𝐺𝑎𝑡𝑒/𝐸𝑙𝑒𝑐𝑡𝑟𝑜𝑙𝑦𝑡𝑒 is relevant on defining 𝐶 𝑡𝑜𝑡 , therefore:

|𝜎| = 𝑉 . 𝐶 𝑡𝑜𝑡 ≈ 𝑉 . 𝐶 𝐺𝑎𝑡𝑒/𝐸𝑙𝑒𝑐𝑡𝑟𝑜𝑙𝑦𝑡𝑒 (38)
Au/MCH interfaces in water can have interfacial capacitances 63 between 1 and 10 uF.cm -2 . However, this greatly depends on the quality of the MCH monolayer, while bare gold in water can reach capacitances 63 up to 35 uF.cm -2 . We fixed 𝑉 at -0.5 V and tested different charges corresponding to this entire range (from 5 to 35 uF.cm -2 ) of capacitance values. These charges were then spread equally among all gold atoms.

Solvent equilibration

For all systems containing the gold surface, we initially perform an energy minimization in three steps: the first step is done by adding harmonic restraints on the gold surface and the MCH monolayer with a force constant of 10.0 kcal.mol -1 .Å -2 in order to keep it in place, followed by a second step where this force constant is decreased to 1.0 kcal.mol -1 .Å -2 and the restraints are removed for the third step.

The heating process is done in two steps: the first one was done in NVT conditions, heating the system from 0 K to 100 K in 50 ps with a harmonic restraint on the gold surface and the MCH molecules with a force constant of 10.0 kcal.mol -1 .Å -2 .

The second step is done in NTP conditions with the harmonic restraints kept in place, heating the system from 100 K to 300 K in 100 ps.

After the system is heated to 300 K, we run a water equilibrium step in NTP conditions for 500 ps in MD simulations on which we completely freeze all the atoms but the free ions and the water ones, not considering bonded energy terms, van der Waals interactions and space electrostatic interactions between any pair of frozen atoms. This is done to decrease the volume of the periodic box on the z direction whereas the surface is on the xy plane, removing any vacuum bubbles that might have formed from the previous steps without excessively deforming the gold surface.

We then perform a second energy minimization step with harmonic restraints of a force constant 10.0 kcal.mol -1 .Å -2 on the gold surface and the monolayer, using as reference to these restraints the final structure obtained from the energy minimization with no restraints previously done.

For the heating process, there is also no longer need to pass through a NTP step, as the periodic box was already optimized from the previous steps, so we heat the system directly from 0 K to 300 K in NVT conditions over 500 ps.

As a side note, production step simulations involving the gold/MCH surface were performed in the NVT ensemble which has been demonstrated as more appropriate for biomolecular systems at solid interface 64 instead of the usual NTP ensemble used in isotropic simulations. All the steps described in this section are schematized in Figure 14.

Analyses of ion distributions

We have performed the counting of Na + , Cl -and phosphate ions using a script previously written by our group 13 that counts the ions and generates a histogram of the heights of the ions in batches of 4 Å, with the origin at the average height of the oxygen atoms of the MCH monolayer up to the Au surface of the periodic box z+1. This is done on each frame of the entire MD trajectory and finally, we manually average this value over the number of analyzed frames. A schematic of this process can be seen on Figure 15. We have also explored the coordination number of the ions in our system in certain simulations to determine the number of water molecules that are present in their first shell of solvation. This number can be compared to experimental results obtained through X-ray diffraction, extended X-ray absorption fine-structure spectroscopy (EXAFS) or neutron diffraction. We have calculated the coordination number from MD trajectories using the radial distribution function ( 𝑔(𝑟) ). By this method the coordination number 𝜂 𝛼𝛽 (𝑟) of atoms of a type 𝛽 inside the coordination shell of an atom 𝛼 extending from 0 to 𝑟 is defined as:

𝜂 𝛼𝛽 (𝑟) = 4𝜋𝑁 𝑉 ∫ 𝑔 𝛼𝛽 (𝑟)𝑟 2 𝑑𝑟 𝑟 0 (39)
on which 𝑁 is the total number of atoms 𝛽 in the system and 𝑉 is the total volume of the simulation cell. We were interested in how the surface and surface charges would influence the solvation of the ions, especially close to the surface, where this number is usually disturbed in the region inside the Debye length of electrodes, that is, the region on which the electrostatic effects of the electrode are at their highest. shown on Figure 1. The closed state presented on the left of Figure 1 is the expected structure when not in presence of AMP, which from now on we will call "off-state", while the open state on the right of Figure 1 is the expected structure in presence of AMP, or the "on-state". The structural formula of AMP can be seen on Figure 2. We have chosen our sequence because it has already been effectively used in an affinity aptasensor by Yu and Lai 2 with promising results, even if it is probably not the most optimal sequence. In this work 2 , the researchers developed an affinity electrochemical biosensor on which the same aptamer sequence was grafted onto a Au/MCH surface by a C6H12-SH function on the 5' extremity and a methylene blue function on 3'. In their system, a conformational switch caused by the presence of AMP changes the distance between the electroactive methylene blue function and the electrode, therefore changing the output electrochemical signal (current in AC voltammetry).

Details on the DNA sequence and conformational switch

The expected conformational switch for the system was designed to answer certain problematics that we will further discuss in this section. First of all, we want the aptamer sequence to be able to easily fold into the recognition hairpin, therefore we decided to put the aptamer sequence on the central DNA strand, which is not directly grafted to the surface, instead only connected to the rest of the system by hybridization (red section on Figure 1). Furthermore, this section of the system containing the aptamer sequence must be readily available to the solvent, therefore a spacer sequence of five adenines has been added between this section and the alkyl thiol section (see Figure 4). For the conformational switch itself, we had two options regarding the dehybridization process: we could either let the central strand de-hybridize on both 3' and 5' extremities or keep one of the extremities in place after AMP binding by adding a very stable sequence of base pairs, such as a GC repeating section. We have chosen the latter option, as it theoretically allows the system to regenerate back to the original state before analysis, after removing the AMP from the medium.

Another important detail is that the structure must be flexible enough to allow the de-hybridization and folding of the hairpin in the presence of AMP. Following this logic, an "induced-fit" -like recognition process is expected, in which the presence of AMP promotes the formation of the recognition loop structure. This type of recognition process has already been reported in literature immobilized aptamers 6 , in larger nucleic acids containing an aptamer 7 and seen in computational studies by molecular dynamics 8 . In all cases, the addition of the target induces a structural change in the aptamer that decreases the free energy of an intermediate structure and directs the process towards the final aptamer-target complex. This effect has been shown by Munzar et al to be capable of overcoming the binding energy of partially hybridized aptamers when those are immobilized to a surface 6 , forming a ternary complex before detaching the aptamer from complementary strand, as seen on Figure 5.

The transduction process in an EGOFET system

One thing to consider about our system is that the conformational changes are likely to take place beyond the Debye length of the electrode if we assume that the measurements take place in electrolytes of ionic force similar to the physiological ionic strength, around 150-160 mM. For this reason, we believe our system is best suited to work in configurations based on dielectric properties, such as field effect transistors (which is our objective) that follow the changes in interfacial capacitance, or optical methods that follow the refracting index.

As DNA is a polyelectrolyte with negative charges in its chain, our hypothesis is that when it undergoes any conformational change, it also changes the overall ionic distribution around it due to the changes in the electrostatic field at the DNA's vicinity.

Therefore, when the DNA is grafted to an electrode and undergoes a major conformational change in the direction normal to the electrode's surface it should directly change the capacitance of the electrode/electrolyte interface. This is further illustrated in Figure 6. As the interactions between the ions and the DNA structure are, overall, of electrostatic nature, we can assume a loose correlation between the sodium ion positioning and the position of the regions of the nucleic acid that are more negatively charged, that is, the phosphate bridges between the nucleotides.

The electrostatic environment around the system can get even more complicated when we consider that the gate electrode in EGOFETs is also electrically charged, therefore also applying an electrostatic force to the ions. As the organic semiconductors commonly used in the channel are of p-type nature 9 , we can assume that most of the time the gate electrode is negatively charged, therefore attracting Na + ions closer to the surface. Following the Gouy-Chapman-Stern model, there are two layers at the vicinity of the interface: there is a compact layer composed by solvent molecules followed by a larger diffuse layer in which the ionic concentration is very high next to the compact layers then decays exponentially until reaching the bulk concentration, as seen in Figure 7. 2. "Off-state" structure construction 2.1 Initial construction of "off-state" DNA structure in solution

The DNA structure was built using the NAB 10 (nucleic acids builder) module of the AmberTools suite 11 . The NAB module works by positioning all bases in order to create the desired nucleic acids helix according to one of the six different types of helix available. The Dickerson conventions of rise and twist are respected in this process 12 .

The backbones constituted by the sugars and phosphates are then added to the structure, which is optimized through an energy minimization process.

The program outputs a PDB file containing the double-stranded B-DNA structure presented in Figure 3 (left and full sequence at the bottom). To obtain the off-state DNA structure, we have manually removed the bases that are not part of it, that is, the bases paired with the single-strand sections of our structure (Figure 8, center and bottom, bases in red) to obtain the off-state sequence (Figure 8, right). Opening this structure on the xleap software while the OL15 force field is loaded automatically corrects any possible missing atom, rendering it ready to be solvated. In order to graft this DNA structure to a Au surface at both extremities, we needed to find an energetically favorable U-shaped conformation on which the extremities are available and pointing at the same direction in space. Our search for this U-shaped conformation is further detailed in the following sections of this chapter, on which we closely followed the evolution of the distance between the 3' and 5' adenines at both extremities of our system (marked with two green arrows on Figure 8).

Initial Molecular Dynamics of the "off-state" DNA structure in implicit and explicit solvent

As a result of the method used to build our DNA sequence, it is initially in an unnaturally straight double-strand conformation (see upper part of Figure 3). In order to properly model the DNA behavior, we need to study its actual conformational space, finding more "natural" (stable) conformations for the sequence in solution and when grafted to a surface. We have performed a series of simulations progressively increasing the complexity of the studies until we considered that the real dynamics of the DNA sequence was being properly modelled.

One aspect of the conformational study we considered was the feasibility of the grafting process of the DNA structure. As the off-state was planned to be grafted by its extremities, one recurrent metric analyzed was the distance between these extremities, as this is directly translated to the distance between the grafting sites on the surface, as well as a direct metric of how far we stray off of the original unnatural straight conformation, as schematized on Figure 9, left. Another complementary metric considered is the radius of gyration of the system (see Figure 9, right). This metric, defined as the radial distance between the center of mass of the structure and a point that as the same moment of inertia as the system if the entire mass were concentrated at this point, is a more direct indication on how folded on itself the DNA structure is, and its value comes from the fact that our system has two very mobile single-strand poly-A tails on both extremities. Therefore, when analyzing the distance between the extremities (Figure 9, left), it is possible that a certain decrease in value might be due to the folding of these two sections without necessarily folding the entire system.

Our first simulations of the off-state system were performed in implicit solvent, modelling the solvent as a continuum media with a defined permittivity or dielectric constant, for a total of 250 ns of simulation while gathering data from the last 200 ns of simulation. This was our first attempt at studying the behavior of the off-state DNA structure, as a more straightforward approach to the conformational study at lower computational costs. We followed the distance between the centers of mass of the adenine bases at both extremities of the system (Figure 10, left) which in the initial structure built by NAB 10 was 145.82 Å (line in red in Figure 10, left). The values obtained did not stray far from the original value over the 200 ns. During this simulation, we have also followed the evolution of the RMSD values of the backbone atoms that fluctuated between 5 Å and 20 Å, remaining between 10 Å and 15 Å with high fluctuations after 100 ns, as shown on Figure 10, right.

These results indicate that we have not properly accessed the conformational space of the DNA structure in implicit solvent. It is known 13,14 that classic models of implicit solvent are not very effective when simulating nucleic acids in general, due to the strong ionic characteristic of these molecules. The results were therefore deemed inconclusive and as the next step in increasing the complexity of the simulations, we proceeded to test the same DNA structure in explicit solvent.

To perform the same simulation on explicit solvent, we constructed the system in a solvent box (TIP3P water molecules) with [Na + ] = 160 mM and Cl -as counter-ions to the final positive charge of the system, due to the negative charges already found in the DNA molecule. We performed 150 ns of simulations on it and analyzed the last 100 ns, this timeframe of analysis being a good compromise between computational costs and exploration. The evolution of the distance between the adenines at the extremities is shown on Figure 11 We can see that around 85 ns, we achieve a minimum value of this distance around 70 Å, however under closer inspection, we can see that this value was achieved mostly by the twisting of the poly-A tails at the extremities instead of the whole system as seen in Figure 12, and therefore was deemed inadequate for the grafting on the surface. These results indicated that the conformational exploration under these conditions was still insufficient. Confronted with this fact, we can:

• decrease the degrees of freedom of our system either by shortening it or passing to a less detailed model, such as coarse grain models,

• Enhanced the sampling with a conformational exploration algorithm, such as REMD or GAMD 15,16 .

Of these, we came to the agreement that insisting on the explicit solvent simulations of the entire system would demand too much of a computational cost, and that the algorithm methods demand an extra carefulness in the assessment of the results obtained 17 .

In our team we recently employed the Gaussian Accelerated Molecular Dynamics (GAMD) technique to enhance the sampling of proteins 18 . However, using this algorithm for nucleic-acids is challenging. For example, GAMD uses dihedral angles energy barriers for smoothing the potential energy and thus facilitate the sampling during molecular dynamics 19,16,17 . To date, only one work reports GAMD for RNA 20 which is the solely example of nucleic acids available on the literature. Coarse

Grain models intrinsically cause a loss of atomic information that we assumed too important for its application in this system. We decided then to decrease the degrees of our system by shortening it.

Considering that our system is mainly composed of two double-strand sections connected by a central de-hybridized section, we can expect this central section of the structure to be able to bend more easily than the double helix sections, and play a strong role on the conformational space of the total structure. This was further confirmed by analyzing the RMSF values of the system referent to the MD simulation in explicit solvent, which is a measurement of the deviation of the coordinates of the particles from a referential average set of coordinates. A peak in the RMSF values per base around this central region was found, as seen on Figure 14, indicating the higher mobility of these nucleotides. Therefore, we came up with the idea of reducing the off-state system to its central section and performing the conformational study only of this region instead, in hopes that this smaller conformational space is a descriptive image of the conformational space of the entire system, obtained at a smaller computational cost.

This process is the main focus of the discussion on the next section.

2.3 Reduction of the "off-state" system and conformational study of the reduced sequence

Our main goal in this study is to develop and find a stable conformation for a reduced DNA structure that we could base ourselves on to generate restraints on the full off-state system in hopes of obtaining a U-shaped structure, assuming that the conformational space of this reduced sequence is a valid and descriptive image of the conformational space of the full off-state system. We have decided to run a conformational study on a smaller DNA sequence composed of the central dehybridized section with four hybridized base pairs on each extremity, as seen on Figure 15 (section in red). We have run six replicas of simulations in explicit solvent of this reduced DNA structure in order to explore its conformational space and analyzed the final 200 ns of each simulation. As the double-strand regions are quite small (four bases each), we have also used harmonic-type energy restraints to keep the base pairs hybridized (5 kcal.mol -1 .Å -2 ). This was done by adding an energetic penalty for all the atoms involved in the hydrogen bond interactions when they are not between 1.80 and 2.04 Å from each other with the energy attaining its maxima at distances at 1.30 and 2.54 Å, for a full penalty of 1.25 kcal/mol. The shape of these potential energy restraints is shown on Figure 16. A similar energetic penalty was also applied for the C1' atoms of the base pairs with optimal distances between 4.90 and 5.00 Å, with the energy attaining its maxima at 4.40 and 5.50 Å. One information obtained from these simulations is that this smaller system has two noticeable main conformations, one open where the main strand remains relatively linear and another where the whole system folds on itself. This can be seen by analyzing the radius of gyration of the DNA structure as seen on Figure 17. The average radii of gyration of these conformations are around 18.7 Å and 14.5 Å, respectively. These two families of conformations, however, are not particularly time dependent and transitions between the families on both directions are observed in the analyzed data, indicating that these two families of conformations have similar energies. In order to further analyze the nature of this reduced DNA structure, we have performed a free energy cartography using Boltzmann statistics as seen on the following equation:

∆𝐺(𝛷 1 , 𝛷 2 ) = -𝑘 𝐵 𝑇 . ln [ 𝑃(𝜌 𝛷 1 ,𝛷 2 ) 𝑃 𝑚𝑎𝑥 (𝜌 𝛷 1 ,𝛷 2 ) ] (1)
For 𝛷 1 and 𝛷 2 , we have tested many different types of structural coordinates, such as distances, angles and dihedrals. When it comes to angles and dihedrals, we found that they failed to give discreet energy minima, instead providing a wide range of values with favorable ∆G values. Instead, inspired on the shape of a telephone handset, we attempted to generate our free energy landscape using two distances between two pairs of nucleotides: D1 being the distance between the phosphorus atoms in the nucleotides T3/G13 and D2 being the distance between the same atoms in the nucleotides C5/C12, as seen on Figure 18. This was done to simulate a hypothetical conformation on the full off-state structure where the extremities are brought as close as possible without overly tensioning the double stranded sections. With these two structural coordinates, we have obtained the free energy landscape which is represented on Figure 19, from which we extracted 4 pairs of values (circles in Figure 19) in local energy minima. These values were to be used as parameters for structural restraints on the full off-state DNA structure further on. Further analyzing the free energy landscape, we can see that a large zone of the free energy landscape is accessible, showing that our reduced system is very flexible. Moreover, we can see that there is a certain linear correlation between D1 and D2 in certain conformation, as there is a region that follows a linear path, indicated with the dotted line. This is to be expected as these two coordinates are not independent and, furthermore, are measured in the same molecule in somewhat the same direction in space when the system in its linear state. However, the regions that interest us the most are the ones on the lower section of the graph, on which D1 has a value lower than the ones found on the aforementioned linear section. Keeping in mind Figure 18, a structure that resembles a telephone handset would have its D1 value decreased With our D1 and D2 values sorted, we were ready to transpose them to the full off-state system.

Molecular Dynamics of "off-state" system with structural restraints

To get back to the structure of the full off-state system (see Figures 8 and15), the values for D1 and D2 found on free energy cartography (Figure 19) were used to define structural restraints (5 kcal.mol -1 .Å -2 ) on the equivalent nucleotides of the offstate structure. After 20 ns of MD simulations with restraints followed by 50 ns of relaxation, we obtained the structures shown on Figure 19. The restraint energy potentials used are of similar nature to the one shown on Figure 16. We have decided to use the structure obtained from the simulation using D1 between 13.50 Å and 15.50 Å and D2 between 12.00 Å and 15.00 Å (Structure A on Figure 19). This choice of structure was done visually and took into account how simple it would be to graft the structure on the surface.

The structure A was lightly retouched in order to make both poly-A tails point at the same direction to be grafted to the surface. This was at first done visually using the Maestro software, followed by a local minimization process encompassing the entire poly-A tail. This last geometry was then grafted through a phosphate linker to the Au/MCH surface of dimensions 9.1 x 8.9 nm on both extremities, forming the off-state bridge system. This was done by removing the hydrogens from the hydroxide group of a MCH molecule as well as the terminal hydroxide group at the end of the DNA sequence, adding a phosphorus atom bonded to two oxygen atoms in between the now available oxygen atoms of the MCH molecule and DNA sequence, binding them to the phosphorus atom to form the phosphate group. The partial atom charges of the new atoms were taken from the phosphates in the DNA sequence and slightly retouched to guarantee that the system maintains an integer as total electric charge.

A water box was added with [Na + ] = 160 mM and Cl -as counter-ions to the final positive charge of the system due to the negative charges already present in the DNA molecule, as shown on Figure 20. Our on-state system, as seen on Figure 21, consists of two different DNA sequences grafted to the Au/MCH surface: one larger DNA structure containing the aptamer hairpin bound to the AMP molecule (Figure 21, left), and a secondary DNA probe (Figure 21, right) that was originally bound to the aptamer sequence while in the off-state. In this section, we will discuss the construction of the on-state system, including the retrieval of the 3D structure of the aptamer hairpin/AMP complex and grafting of both DNA structures to the surface.

Elucidation of the AMP-aptamer complex structure

In order to build the on-state DNA structure bound to AMP (Figure 21, left), we first elucidated the structure of the aptamer-AMP complex. To date, there is no experimental 3D structure of this aptamer-AMP complex. Song et al first published three sequences of anti-AMP aptamers which were grafted on magnetic beads and used in a fluorescence-colorimetric strategy 5 ; that is the only publication found where anti-AMP sequences were developed. Out of these three sequences, we have chosen the same sequence that was used by Yu 2 in an aptasensor that was also immobilized to a gold surface, although using a different type of transduction method.

This elucidation was done using a computational pipeline method previously published by Jeddi et al 21 on which the sequence is first analyzed using the mfold 22 platform in order to obtain the secondary structure of the hairpin. When we did so, the results were in consensus with the ones proposed by Yu and Song (Figure 22a)). Once we had this hairpin structure, we used it to generate 3D structures of it as RNA on the 3dRNAweb 23 platform. We took the five best scorings obtained from 3dRNAweb, manually converted it back to DNA and ran relaxation simulations for 50 ns, obtaining the structures shown on Figure 22b). We have engaged in molecular docking computations using the AutoDock 24 software, defining the entire hairpin structure as interacting site therefore rendering the docking process unbiased when it comes to choosing the binding cavity. Nevertheless, one common binding cavity is found for the 5 aptamer structures, as shown on Figure 22c) even though this ligand could interact with any other region of the DNA structure. When Song et al first proposed the three anti-AMP sequences, the authors conjectured that a common sequence of GGTT in the vicinity of a GC base-pair could be the interacting site of the AMP. The GGTT structural loop is not exactly at the cavity but instead adjacent to it, with the cavity around the GC base-pair. While the interacting cavity on the hairpin is roughly similar for all aptamer structures, recurrently formed by the bases fromT6 to G9 on one side and the bases from G15 to G18 on the other side.

This cavity is shown in more detail on Figure 23. The interacting modes found show the AMP interacting with this cavity in different orientations on each structure. The ΔG values obtained by the molecular docking of the five structures are shown on Table 1. To test the stability of each structure, supplemental molecular dynamics simulations for 100 ns in water boxes were realized for the five complexes. Out of all the structures tested, only Structures 2 and 4 kept the AMP in place over 100 ns, both are represented on Figure 24. We have used the MM-GBSA method to compute the binding energy of both complexes and found two favorable energies of ΔG2 = -41.4 ± 6.2 kcal.mol -1 and ΔG4 =-35.6 ± 4.0 kcal.mol -1 for structures 2 and 4, respectively. On both structures, the free energies are of very similar values that actually have overlapping uncertainties.

In order to define the most structure for our further studies, we have further analyzed both of them visually. Structure 2 was chosen due an additional hydrogen bond formed between the aptamer and AMP and the amide bond in AMP adopts a favorable trans conformation, which is not the case in Structure 4, as also shown on Figure 24. For all the aforementioned reasons, Structure 2 was selected for the construction of the onstate system.

Preparation and grafting of the "on-state" structure

The aptamer/AMP complex structure 2 found in the previous section was manually linked to the rest of the on-state main DNA structure as shown on Figure 25. As this linkage was to be done at the very center of the hairpin, the eventual adaptations on the last base's (G19) conformation were kept to a minimum. We have also added an extra T base at the 5' terminus of the system for further stability. The rest of the DNA structure was built using NAB 10 in the same method as previously discussed in section 2.1. Again, this method generates an unnaturally straight B-DNA double helix structure, so a relaxation step of 50 ns in explicit solvent was performed on this new system. For informative purposes, the RMSD plot of this step is shown on In this relaxation step, we could see a structuration of the hairpin towards the same direction as the double helix, as seen on Figure 27, leading the complete structure to a linear conformation. This seems to be further stabilized by a π-stacking interaction between the extra T base we added at 5' and the G base at the 3' extremity of the other DNA strand that will be grafted to the surface (see Figure 27c)). This was confirmed by the quantification of the non-binding interaction energy (summation of the electrostatic and van der Waals interactions, see chapter 2, section 1.2) between the two residues with respect to the distance between the C1' atoms of each residue (Figure 27d)). and the AMP molecule (green) at t = 0 ns of the relaxation step. b) Snapshot of the same structure at t = 50 ns. c) Detail on the b) structure putting in evidence the π-stacking formed between the 5' terminus of the hairpin and the 3' terminus of the hybridized strand to be grafted to the surface (see Figure 25). d) Measured interaction energy between the two bases plotted against the distance between the C1' atoms of each base.

The final structure was grafted onto the Au/MCH surface of dimensions 9.1 x 8.9 nm 2 by a phosphate linker in the same process described at section 2.4. The secondary DNA leg (Figure 21, right) which is de-hybridized from the bridge (5'-AACCCCAACTAAAAAAAA-3') was also generated by NAB 10 and relaxed by MD simulations in explicit solvent for 50 ns. The RMSD values of this relaxation process is

shown on Figure 28. The final structure obtained in this process was then attached to the same Au/MCH surface by another phosphate linker by its 3' terminus. The final structure constituted by the surface with these two attachments constitutes the on-state system and a solvent box was added with [Na + ] = 160 mM and Cl -as counter ions (due to the negative charges already present in the DNA molecule)

to the final positive net charge (Figure 29). For the sake of clarity, water molecules are not displayed. Figure 31 shows the initial and final frames of one replica of both off-and onstates. No collapse of the DNA structure onto the surface was observed in any simulation. This is in agreement with the experimental results available in literature, that say that MCH monolayers tend to make nearby DNA strands take an upright conformation 25 , this upright conformation was seen especially on the 3-4 bases closest to the surface after a certain simulation time. In addition, we found that the DNA base pairs remain hybridized for the entire analyzed timeframe, emphasizing that the chosen sequences are stable enough even after the grafting process. The ampicillin remains bound inside the binding cavity in the on-state system demonstrating the high affinity of the aptameric sequence and the sturdiness of the hairpin. Finally, for the on-state simulation, there is no re-hybridization of the DNA strands reforming the off-state.

These results reaffirm the stability of our system on both states, which again, theoretically suggests a more sensitive biosensor by reducing the noise in the output signals.

Studies on the impact of surface grafting

To evaluate the influence of the surface immobilization, we compared the MD simulations performed for the systems previously presented to those performed in solution. Therefore, the on-state DNA structure described on Figure 25 was further investigated in solution for comparison with the complete on-state system on surface for 100 ns. The RMSF values were computed for all MD trajectories in order to investigate the variation of flexibility between the systems in solution and grafted onto the MCH gold surface. For the on-state, we have found a slight decrease of relative flexibility near the grafting site, as evidenced in the RMSF values per residue shown on Figure 32. For the off-state, the RMSF values decreased when compared to the values found in solution, as illustrated on Figure 33. This is accentuated on the extremities which are directly linked onto the surface, and there is also a sharp decrease of the RMSF values for the bases localized at the central de-hybridized DNA section. This

shows that the effects of the grafting are not only felt by the nucleic acid's residues close to the grafting site but also by the whole DNA structure, demonstrating the need for paying extra attention to the conformation chosen for the grafting, as well as the importance of choosing an optimized sequence. There is also a sharp decrease in the mobility of the central de-hybridized section, reaffirming the importance of the conformational study we performed in the reduced system representing this section (described in section 2.3 of this chapter).

Analysis of ionic movement and comparison of ion profiles between states

In the previous section, we have shown how feasible our DNA system is for the grafting and its stability in both off-and on-states. From a mechanical point of view, our system is successful. However, to go further in the discussion, we want to translate our results in Molecular Dynamics towards the actual output signal of the EGOFET-based biosensor. As a reminder, a critical point of the transduction process in an EGOFET biosensor is the change of the capacitance at the interface between the gate electrode and the electrolyte. This capacitance is directly influenced by the organization of the ions in the solution at the interface of the electrode, according to the Gouy-Chapman-Stern model (see Figure 7). To simulate an ionic force similar to the usual experimental conditions, we have performed all simulations with the addition of Na + ions up to the concentration of 160 mM and Cl -ions up to the system neutrality, due to the negative charges already present in the DNA molecule. This ionic concentration is commonly described in literature for molecular modelling as it is similar to the ionic strength found in physiological media.

As a way of connecting the proposed conformational structures to the interfacial capacitance and therefore the output current signal of an EGOFET biosensor, we analyzed the ionic distribution over the last 100 ns of each replica by counting the ions over slabs of width 4 Å, parallel to the surface. This width has been previously found by our team as the most reliable one for this type of analysis of ionic distribution 26 . To obtain an average atom count, the total ionic count was divided by the number of analyzed frames (20000 frames for 100 ns).

We analyzed the Na + concentration profiles in both off-and on-states and observed a noticeable difference. The regions on the concentration profile of Na + with a noticeable increase in value extend up to around 8 nm and 10 nm from the surface for the off-and the on-state respectively, with the off-state concentrating it in a single distance range while the on-state presents it in two subsequent distance ranges. There is a direct link between the concentration profiles of Na + and the phosphates groups in both states, indicating that the ionic distribution change is directly connected to the architecture conformation switch upon the AMP binding and the Na + ions position is mostly governed by the electrostatic interactions between the ions and the DNA structure, which is, as all nucleic acids, a polyelectrolyte. The secondary peak in the phosphate profile is also related to the position of the aptamer hairpin, where the phosphates are tightly packed. All these profiles are available in Figure 34. To properly evaluate the effect of the DNA structure on the ionic profile, we have also performed MD simulations on a Au/MCH surface without any grafted DNA for comparison. Figure 35 shows the Na + and Cl -distributions of such system. Both Na + and Cl -concentration profiles have the same feature and are homogeneous above a distance to the surface of approximately 2.5 nm. The Cl -concentration profile for our systems containing our DNA structure shows a very similar behavior (with just a slight decrease in the regions containing the DNA molecules) to the surface without DNA, meanwhile the DNA effectively increases the Na + concentration in the vicinity of the surface for both states (see Fig. 34 c)). In crystals, the coordination number is the number of ions that immediately surround an ion of the opposite charge within a crystal lattice. In water solution, this metric computes the number of atoms from water molecules directly surrounding the ion, as illustrated in Figure 36. Since ions are the key elements of the detection of our biosensor, we investigated if their intrinsic behavior are affected by the surface and the DNA. To do this, we computed the coordination numbers of Na + and Cl -ions. Therefore, we analyzed the coordination number of the ions in both off-and on-states, with the results of this study shown on Table 2 compared to ab initio results described in literature 27 .

The equation used to calculate the coordination number of an atom 𝛽 in the coordination shell a central atom 𝛼 is the following (this is further discussed in Chapter 2 section 6.4):

𝜂 𝛼𝛽 (𝑟) = 4𝜋𝑁 𝑉 ∫ 𝑔 𝛼𝛽 (𝑟)𝑟 2 𝑑𝑟 𝑟 0 (2)
on which 𝑔 𝛼𝛽 (𝑟) is the radial distribution function of 𝛽 with respect to 𝛼, 𝑁 is the total number of atoms 𝛽 in the system and 𝑉 is the total volume of the simulation cell.

To discriminate the ionic environment in the vicinity of the device and keeping in mind the Gouy-Chapman-Stern model, we have defined our diffuse layer between 0 and 8 nm due to the influence of the negative charges in the DNA molecule through their phosphate sections and bulk as the solution that is outside the range of the DNA structure, that is, over 8 nm away from the surface. We have also used the same cutoff values as the ab initio study 27 for the integration on the radial distribution function (𝑟, as the integration upper boundary): 3.2 Ǻ for Na-O, 3.7 Ǻ for Na-H, 3.8 Ǻ for Cl-O and 2.9 Ǻ for Cl-H 27 . When comparing the obtained results from our system, it shows similar values on both diffuse and bulk regions than the results from the literature. On top of that, the results are also in agreement with experimental results for the Na-O number, being between 4 and 8 28 . These results validate our choice of TIP3P as an effective water model.

Even if the coordination numbers are slightly lower in the diffuse layer than in bulk, this difference is not significant. Other defined distances of diffuse layer much closer to the surface were tested and provided similar results. Therefore, we conclude that our system does not induce a significant perturbation of ionic solvation. We also conclude that these slight decrease in the coordination number comes exclusively from the DNA structure, as if we consider the theoretical Debye length of 7.71 Ǻ for our solution and the length of the SAM of around 7.6 Ǻ, the effective distance inside the Debye length is only around 0.1 Ǻ above the SAM. This, added to the fact that the water and ions rarely penetrate the SAM during the simulations, leads to the conclusion that all these measurements are done beyond the Debye length. Another important conclusion from these results is the lack of any major differences in the ionic solvation between off-and on-states.

The observation of no major differences of sodium coordination number in the diffuse and bulk layers is particularly interesting. Indeed, the popular belief is that sodium ions coordinate to the phosphate moieties of DNA 29 . In reality, this basic idea came to the fact that monovalent ions are difficult to distinguish from water molecules in crystallographic studies 26 . In 2014, Richard Lavery and collaborators made long trajectory of molecular dynamics (around µs scale) and elucidate a more subtle role of monovalent cation around the DNA helix 30 which are, of course, around the phosphate moieties but without a definite coordination. Indeed, the nucleic acid acts more as an electrostatic field and attracts ions without desolvating them. In our work, we observe that Na + attracted by DNA (see previous chapter) but the fact that the Na + coordination number does not decrease by, at least, one unit highlight this subtle role of electrostatic field of DNA.

To sum up all the results shown in this section, the proposed conformational switch would effectively change the surface capacitance (see Figure 6); therefore, the system is a valid biorecognition component for a FET-type AMP biosensor, although we have found no significant changes in the ionic solvation near the surface.

Testing the effects of sequence length on the ionic distribution

In the previous section, we showed that our off-and on-states have effectively different effects on the ionic distribution next to the surface, especially on the Na + ion profile, that seems to follow the position of the phosphates on the DNA sequence, although not by strict coordination of these negatively charged centers. To further optimize the system, we have tested the effects of sequence length on the ionic distribution over the surface. As the conformational study described on section 2.3 has a considerably high computational cost, we decided to avoid extensively changing our system so we could still use the results from it. So, the approach we took for this test was to remove the extra T base we added next to our aptamer sequence at the 5' terminus on section 2.1 of this chapter and follow up by removing another T base, this time from the aptamer sequence itself, forming the mutant sequences shown on Figure 37. We have run two replicas of MD simulations of 100 ns on these two new sequences on the on-state only, as we deemed that due to the structural restraints of the off-state at the changing section of the DNA structure, simulating these new sequences at the off-state would be of no interest. Figure 38 shows the Na + ion profile for these two new sequences in comparison to the original sequence. When comparing the Na + profiles between the original base sequence and the mutant sequences, one noticeable change is the diminution and approaching of the secondary peak towards the main peak. As we had previously established in section 4.3, this peak is directly related to the position of the aptamer hairpin, so we decided to further investigate the structure of this hairpin during the simulations. What we found was that the shortening of the sequence rendered the previously discussed π-stacking interaction between this strand's 5' terminus and the support hybridized strand's 3' terminus less effective, as we can see in the snapshot presented in Figure 39. With this interaction weakened, the aptamer hairpin does not stabilize at an upright position compared to the rest of the DNA structure. This effect is further amplified by the fact that the rest of the structure doesn't necessarily stay upright during the entire simulation, bringing the average height of the aptamer hairpin even closer to the surface. Even though this lowered average height of the aptamer hairpin brings the second peak in the Na + ion profile closer to the first peak, we are not certain that this is the only cause for this result. It is also important to remember that the main interaction between the Na + ions and the DNA structure is electrostatic, and that removing bases from the DNA structure effectively removes charges from it, therefore lowering its overall electrostatic force applied on the Na + ions. We do believe, however, that the relative position of the DNA structure is the main defining aspect of the shape of the Na + ion profile.

This result also highlights the importance of an upright DNA structure at the onstate for an effective difference in ionic distribution between the states. This must be taken into consideration for all future design processes of this type of aptasensor. We propose that the most efficient conformational switch for this aptasensor would be that the off-state conformation remains as close as possible to the electrode while the onstate conformation extends as upright as possible, distancing itself from the electrode, as represented on Figure 40. Another conclusion of this study was that our initial sequence gave better overall results, therefore we continued using it in all further studies.

Evaluation of the hybridization energy of the aptameric sequence

The transition between the off-state and the on-state of our system involves the de-hybridization of the aptamer sequence from the rest of the DNA structure and the formation of the recognition hairpin. Although a study of this entire process would go way beyond the scope of this thesis, we have evaluated the hybridization free energy of the aptamer sequence by Umbrella Sampling (US) to give a first insight on the energy values over the transition between states.

As reaction coordinate, we have chosen the collective distance between the atoms C1' of the nucleotides of the aptamer and the C1' atoms of the complementary strand that is hybridized to it at the off-state, as seen on Figure 41.

Validation of the umbrella sampling method

At the end of the umbrella sampling calculation, like for any simulation procedure, we must check several items in order to ensure the validity of our results.

First of all, according to the definition of our reaction coordinate, we visualize the simulations made with the collective distance around 40 Å. We observe that the strands are de-hybridized as illustrated on Figure 42 indicating this way that the reaction coordinate setting was correct. In this chapter, we focused on the off-state system, that is, the system grafted to the MCH gold surface. Similar results were found for the system in solution. One of the most important parameters of a successful umbrella sampling computation is that the reaction coordinate overlaps between adjacent windows, so that the reaction coordinate is entirely explored. This requirement is extracted by analyzing the results with the Histogram_plot software 31 which provides us Figure 43, a plot of the population of each window as a function of the reaction coordinate, where we observe no gap between windows and an effective overlap. We also observe a Gaussian distribution of windows which is in accordance with the harmonic potential implemented in the umbrella sampling technique. The good overlap of windows is also observed by looking at the phase space covering. This one is obtained through the help of the Discrete_plot software 32 and display on Figure 44. Here, the complete covering indicates that there is no unexplored area of the reaction coordinate, as expected. as in solution. This low hybridization free energy variation when the nucleic acid is grafted on the surface could be an asset for our future device.

However, the curve profile is not exactly similar for the system when it is grafted on the MCH gold surface as the free in solution. In fact, the dissociation of DNA in solution present three energetic phases along the reaction coordinate:

• A steep increase when the collective distance is raised from 8.8 to 13.3 Å.

• A flat part between -8 kcal/mol and -6 kcal/mol for reaction coordinates from 13.3 to 31.8 Å.

• Another steeper increase curve from 31.8 to 40 Å in order to reach the complete dissociation.

In contrast, when the DNA sequence is immobilized on the surface, we observe a steadily rise of the free energy when the reaction coordinate increases.

To elucidate this discrepancy, we observe some structural features of the DNA in solution with various reaction coordinate values. Some typical ones are displayed on Figure 46. In fact, when the DNA separation distance rise from 13.2 to 32 Å, the DNA segment find ways to accommodate the restraint and undergoes various reorganizations. This fact highlights the high plasticity of DNA in solution 34 . Indeed, contrary to the popular belief of a canonical DNA helical structure, this biomolecule is flexible and proposes many structural organizations 35 . Other observations can be made from Figure 46. The dissociation starts from the 5' extremity of the aptameric sequence which correspond to the part near the surface when the system is grafted on the MCH gold surface. This fact is a good prefiguration of our aptasensor and will be hereafter confirmed on the umbrella sampling simulations of the grafted DNA. The other interesting fact is that the aptameric sequence, when starts to be dissociated, does not remain as a single strand but rather undergoes a folding. This prefigures the formation of the hairpin structure of the aptamer structure even if this last structure is far from being identical to the functional aptamer. When the DNA is immobilized on the MCH gold surface its dissociation mechanism is slightly different. Figure 47 illustrates typical structures issued from US simulations and with various values for the reaction coordinate. As witnessed by the PMF curve of Figure 45, the loose of interaction is regular as the reaction coordinate increases. Like for the DNA in solution the dissociation starts from the 5' extremity of the aptameric sequence (near the surface). However, conversely to the DNA in solution which steadily reorganize its structure to compensate the increase of reaction coordinate, the aptameric sequence slips along the DNA fragment immobilized on the surface. At the end of the simulation, when the aptameric sequence is dissociated, we do not observe any prefiguration of a hairpin structure as for the DNA alone in solution.

This folding process is now under investigation through Gaussian Accelerated Molecular Dynamics simulations (GAMD) and proposed as prospect of this work.

Although the results shown in this thesis are very preliminary, we hope that they instigate further investigations on the dissociation process.

Evaluating the impact of electrically charged surfaces on the ionic profiles

In previous sections, we tested our system on both off-and on-states by MD simulations and we concluded that our system is stable on both states. We have also shown that we obtain different ionic distributions next to the surface depending on the current state of the system. We also optimized the sequence length for maximizing the variation of ionic distribution and we performed a preliminary evaluation of the dehybridization energy barrier that the system needs to overcome to transition between the off-and on-states by Umbrella Sampling.

In this section, we want to study evaluate the effects of an electrically charged surface on the system, both on the DNA structural variation and in the final ionic distribution. To do this task, we first need to translate the experimental electrical conditions to the atomic scale as used in our models.

7.1 Estimating the electric charge at the surface on a hypothetical

EGOFET system

For our MD simulations, we have chosen to mimic the effects of an applied voltage at the electrode by adding the negative charges directly on the gold atoms.

The process of choosing the right electric charges to add to each atom was discussed in section 6.2 of Chapter 2. As a quick recapitulation, we take into account the system acting as a capacitor and we consider the charge at the surface of the capacitor as the product of the applied voltage V times the capacitance C (Q = C.V).

The value of V can be chosen experimentally, but it is reasonable to assume the potential difference applied during the analysis would be around -0.5 V (reminding that EGOFETs operate under gate voltages under 1 V, and that most organic semiconductors are p-type, i.e., operate at negative gate voltages), therefore we will use this value of |V| = 0.5 V.

We have considered that a gold surface covered by a MCH SAM can have capacitance values between 5 and 35 uF.cm -2 according to literature 36 . On a surface of dimensions 9.1 x 8.9 nm 2 , those values of capacitance and voltage translate to charges Q with a range between 12.65 and 88.58 e in total. As we needed integer values for our total system charge, we have decided to run tests adding charges of 13, 25, 38, 50, 63, 75 and 88 e. These values would be spread over all the gold atoms in the systems (2220 atoms in total). Overall, as these values are defined as C.V rounded to an integer, they account for the large range of capacitance. The correlation of these values and the capacitances can be found on Table 3. These values were added to both on-and off-state systems, adding up to 14 different systems (seven values of charge and two systems for each value).

Effect of surface electrical charging onto DNA conformations

We have manually added the charges to Au atoms found on Table 3 in the structures we had already built in the previous sections and we ran 150 ns of simulations in explicit water and analyzed the last 100 ns. We have performed the simulations in at least five replicas for each value of electric charge added to the surface on both on-and off-states for a total of at least 10 replicas per value of charge.

One issue found from the simulations is that due to the periodic boundary conditions applied, part of the ions gets trapped by the electrostatic force applied by the surface at the image periodic box z+1, therefore the baselines effectively change as we change the electric charge at the surface, as more charge implies more ions trapped. This can be seen on Figure 48, which is a screenshot of one of the replicas with total surface charge 88 e, where this effect is the most noticeable. We attempted to minimize this effect by enlarging the periodic box in the z-axis, therefore rendering the z+1 surface farther from our region of interest, with our new boxes attaining around 26 nm in size over the z-axis. Otherwise, our system remained stable even under the harshest electrostatic conditions, no de-hybridization was detected nor collapse of the DNA structure over the surface, as seen on Figure 49, which presents screenshots of the systems after 100 ns of simulation at the highest tested surface charge of 88 e. In this figure we can see that all the remarks about the surface stability are still valid: our off-state does not de-hybridize, our DNA structures do not collapse over the surface, the aptamer hairpin remains well-structured and the AMP molecule does not leave the interaction cavity.

All these facts emphasize the quality of our chosen DNA system and its robustness.

Ionic profiles of both states under different electric charges

Analyzing the Na + profiles of the simulations, we came to notice that there was a vast section of our box where the Na + concentration remains constant, between 12 nm and 20 nm, as shown on Figure 50, exemplified by the profiles of both off-and onstates with a surface charge of 50 e. In Figure 50, we can also see the region where the ions are trapped by the surface z+1 at around 25 nm, with a huge peak in the profile.

This artefact comes from the phenomenon shown on Figure 49, and to take it into account in our results we needed to perform a normalization of our profiles in different charges. This artefact is also the reason why for the following analyses between the different applied electric charges, we will only analyze the ion profiles up to 12 nm of distance from the surface. We have taken advantage of the region of constant Na + concentration and used it to define a baseline when comparing the profiles of different charges, as the quantity of ions trapped at the end of our box is dependent on the electric charge on the surface.

The overlaying Na + profiles of all the electric charges tested can be seen on Figure 51. The peaks originally found with an electric neutral surface are still present in all profiles, with a major difference being a third peak appearing next to the MCH monolayer due to the increasing electrostatic attraction from the Au molecules. nm for that specific profile. All profiles were done from at least 5 replicas.

To further confirm that the surface charges equally influence both the off-state and the on-state ionic profiles and therefore does not interfere with the difference in ionic distribution between the two profiles, we have also come up with the simple metric 𝛥𝐼𝑜𝑛𝑠 of change in the ionic profile by the following equation:

𝛥𝐼𝑜𝑛𝑠 = 𝑎𝑡𝑜𝑚 𝑐𝑜𝑢𝑛𝑡(𝑜𝑛𝑠𝑡𝑎𝑡𝑒) -𝑎𝑡𝑜𝑚 𝑐𝑜𝑢𝑛𝑡 (𝑜𝑓𝑓𝑠𝑡𝑎𝑡𝑒) (3)
In Figure 52, we present 𝛥𝐼𝑜𝑛𝑠 for all electric charges we tested in our surfaces.

As it can be seen on this figure, this new metric that directly measures the difference between the average atom counts in the on-state and off-state present the same profile for all charges, with an initial section on which the off-state structures show a higher ionic concentration closer to the surface followed by section on which the on-state shows a higher ionic concentration and finally arriving on the "bulk" concentration that is roughly the same between states, therefore showing that even at the highest surface charges our system is still able to show different ionic profiles between on-and offstate. With these results, not only we can conclude that the system is stable under electric fields expected in an EGOFET device, we can also affirm that our DNA structures have a more impactful effect on the ionic distribution than the electric field induced by the electric charges at the gold surface. These electric charges are not enough to cause our systems to display similar ionic distribution profiles next to the system, which ensures the change in interfacial capacitance when the system changes states and therefore the transduction process of our EGOFET biosensor.

We have also checked the hydration coordination number of the ions in solution when our surface is charged with 75 e in off-state, giving the results in Table 4, comparing it to ab initio results in literature 27 , and diffuse layer defined between 0 and 8 nm from the surface and bulk the region after this layer ( > 8 nm away from the surface). When comparing these results with the ones shown on Table 2, we can see that the charges on the surface had no effect in the ionic solvation. Again, we have tested different distances for the diffuse layer with no significant change. This fact induces that the presence of an electrostatic field does not impact the ionic water coordination behavior. Besides, this study of various external electrostatic fields confirms, in a way, the local role of DNA as electrostatic field since there are no discrepancies of water coordination for the Na + in the bulk or diffuse layers.

General Conclusions and Perspectives

We have shown our results on the construction of our model biorecognition structure in both off-and on-states, representing the system in the absence and presence of ampicillin (AMP), respectively. We have effectively found stable conformations for both states, including the elucidation of the complex AMP-aptamer structure. We have also shown the results of the computational studies we performed in the system:

• The system is robust enough for the time measurements and showed no sign of denaturation over the simulation time periods, including de-hybridization of the DNA double helixes nor deformation of the aptamer hairpin. The AMP molecule also stays in place inside the hairpin during the simulations.

• We have analyzed the ionic profiles at the interface in both off-and on-states and we have concluded that they are effectively distinct and that would indicate that these states have distinct interfacial capacitance, therefore, generating different output signals of an EGOFET system whether in presence or absence of AMP. We have also followed the ionic solvation in our systems, although these results showed no difference between diffuse layers and bulk.

• We have tested the effects of mutations by deletion in the ionic distribution and we have performed a small optimization on sequence length. We have also analyzed by Umbrella Sampling the de-hybridization free energy of the aptamer sequence from the rest of the off-state structure and compared it to the de-hybridization of the system in solution, showing the changes in the PMF between the grafted and in solution systems, which indicates a slight change in the de-hybridization mechanism caused by the grafting.

• We have evaluated the effects of surface charge in our DNA structure and ionic distributions: the DNA structures are stable even under harsh electrostatic environments, and apart from a sharp increase in the Na + concentration close to the surface, the ionic distributions are still distinct between off-and on-states.

Graphs of the Na + ionic profile for each value of added electric charge on the surfaceon-state:

Graphs of the Na + ionic profile for each value of added electric charge on the surfaceoff-state:

Résumé en français:

INTRODUCTION

Aptamers are structured DNA or RNA sequences selected in vitro to bind specifically to a non-nucleic acid target. 1 Current SELEX (systematic evolutiono fl i g a n d sb ye x p o n e n t i a l enrichment) methodologies allow us to quickly identify specific aptamer sequences against a broad variety of targets, for example, from whole-cell to low-molecular-weight molecules. Aptamers have proven their relevance as specific bioreceptors and used to date in many highly sensitive analytical methodologies. 2,3 Many of them, named aptasensors, involve the anchoring of one of the aptamer ends to a surface acting as the transducer element. 4 Among the varieties of aptasensors, those that do not require labeling of the distal end of the sequence (i.e., by a biomolecular probe, a protein, a nanoparticle, etc.) to assay the target constitute the most interesting ones according to their ease of use and their often better analytical characteristics. This is especially true for lowmolecular-weight targets that can show dramatic affinity decay upon even minor aptamer labeling. 5 To avoid labeling, powerful transduction strategies have been proposed, those based on the monitoring of the charge distribution change close to the surface upon target binding being particularly promising. Nucleic acid sequences are polyelectrolytes and their conformational change associated with the target binding can lead to charges in reorganization close to the transducing surface. Current analytical techniques such as surface plasmon spectroscopies and field-effect transistors can detect small ionic distribution changes of an electrolyte. The intensity of the recorded signal depends on the refractive index and on the capacitance of the thin solution layer close to the surface, respectively. These techniques have already shown their relevance in the field of aptasensors. 6,7 However, the inherent lack of affinity of aptamers selected against low-molecularweight targets, with affinity constants usually in the micromolar range, requires innovative amplification strategies to reach lower limits of detection.

Recent advances in designing functional DNA architectures 7,8 allow us to consider a new class of highly sensitive aptasensors that switch their entire conformation even if only few bases are involved in a target binding process. In most of the cases, the aptamer sequence is introduced inside a host matrix (scaffold) participating in its overall three-dimensional (3D) structure. 9 By properly designing the nucleic acid scaffold, that is, when anchoring aptamer functions to the scaffold are the same as those interacting with the target, scaffold structure change may be induced by the target binding. The equilibrium constant of the competitive exchange reaction involving the aptamer sequence, the scaffold, and the target must thus be carefully tuned in order to maintain both a scaffold stability at ambient temperature and a relevant apparent aptamer/target affinity constant.

To achieve the design of such structures, information about the structural mechanism of the biomolecular complex formation is required. However, the precise aptamer/target complex structure is known only for a few systems. In these well-described cases, the aptamer is already folded to bind the target constituting a relevant starting point for the overall structure determination. 10,11 However, in the most general case, "induced-fitm e c h a n i s m " 12-14 is involved and the structure cannot be single-handedly deduced from classical tools such as mfold alone. Besides, the vast majority of biosensor devices require immobilization of a biomolecule on a surface, which involves surface/biomolecule interactions and thus structural modifications. Despite its wide importance, experimental investigations of such structural modifications remain elusive so that molecular simulations, 15 in this context, represent an alternative to understand the underlying physical processes which occur at the biomolecule/surface interface. 16 The available literature on molecular modeling studies of sensing architectures including an immobilized bioreceptor is scarce. Some recent articles deal with peptide-based systems (i.e., ochratoxin, 17 carboxyhemoglobin, 18 and dopamine 19 ), but very few engage a nucleic acid sequence as a receptor. Puangmali and collaborators have published a sequence of articles in which they deciphered the mechanisms of naked-eye colorimetric aptasensors, first for the detection of ochratoxin A, 20 then on cysteamine-capped 21 and citrate-capped 22 gold nanoparticles for 8-oxo-dG detection, using molecular dynamics (MD) simulations, with the latter two also showing experimental validation. Drew and Walsh also recently published the study on the effects of temperature and salt concentration on the adsorption of a DNA hairpin onto a gold surface. 23 The modeling of the most encountered aptasensors format, where the nucleic acid sequence is immobilized onto a solid, was only discussed, to our knowledge by Ruan et al. They analyzed the effects of adding a T6 linker to an antiflavin RNA-based aptamer and its anchorage to a 1-mercapto-6hexanol-functionalized gold surface on the biorecognition process. 24 They concluded that the immobilized aptamer undergoes a structural reorganization, which may cause a loss in affinity, and that this restructuring is due to an anisotropic ion distribution around the aptamer. This loss of affinity is dominant when immobilizing bioreceptors to a surface, which is well known experimentally. 25 Our aim is here to take advantage of this ionic distribution variation, at the vicinity of an aptamer grafted on a gold surface, to create an efficient biosensor.

In this work, we propose a DNA-based biorecognition structure immobilized on a gold surface, taking advantage of a conformational switch to promote a change in the ionic concentration at the vicinity of the gold surface, which can be a label-free method of transduction for a future EGOFET-based biosensor. The recognition process is mediated by the presence of an aptamer sequence available for interaction with the target molecule. As a proof of concept, we have chosen ampicillin (AMP) as our molecule of interest, which is environmentally important as it plays an important role in bacterial antibiotic resistance. 26,27 As aptamers can be readily developed for a large group of molecules, we believe that the results obtained with this proof-of-concept would be straightforwardly adapted to other similar analytes.

COMPUTATIONAL DETAILS

2.1. Aptasensor Description. Our system has been designed while keeping in mind its experimental feasibility. Alkyl-thiol self-assembly on gold 28,29 is one of the most popular methods to immobilize a DNA structure in an aptasensor. Grafting can be easily done using 5′ or 3′ alkylthiol-modified DNA sequences. Moreover, saturation of the gold surface with a short alkyl-thiol such as mercaptohexanol (MCH) is commonly realized in order to displace occasional impurities and avoid the collapse of biomolecules onto the uncovered surface. As a consequence, we investigated the structure depicted in Figure 1.

The supramolecular DNA structure is composed of two sequences immobilized on a gold surface through the Au-S bond and bridged through a third strand that contains the anti-AMP aptamer sequence at one of its hybridized extremities (off-state, see Figure 1 left). The presence of AMP in the media is expected to trigger the dehybridization of the aptamer sequence in order to fold with the ligand, leading to another supramolecular structure (on-state, see Figure 1 right).

2.2. In Silico Constructions. DNA systems were built with the nucleic acid builder module of Amber software. 30 Manual edits were performed to generate the sequences displayed in Scheme 1. The initial aptamer hairpin structures were obtained with the help of the mfold 31 platform and 3D RNA web server. 32 The gold surface model used in this work is the Au{111} functionalized with a 6-mercapto-1-hexanol selfassembled monolayer as already defined and validated in a previous work on aptasensor simulations. 24,33 This surface containing 2220 atoms over two atomic layers is of size 9.1 × 8.9 nm 2 , which is large enough to encompass the off-or onstate. Each periodic cell contained one of these surfaces on the xy plane, with a saturation coverage c(2 × 2) for the SAM, totaling 350 molecules of MCH per box, giving us a MCH grafting density of around 4.32 molecules per mm 2 . We have grafted one copy of our proposed DNA structure per cell by grafting the sequence 5′-AACCCCAACTAAAAAAAA-3′ by its 3′ extremity and the sequence 5′ -AAAAAAGCGCGCGCGCGCG-3′ by its 5′ extremity with a distance of 7.2 nm between each grafting site. Due to the labile nature of the Au-S bond, 34 it is hard to determine the appropriate distance between the grafting points; instead, we chose this distance based on a free energy approach we discuss The Journal of Physical Chemistry B further in the Results and Discussion section and practicality in the system construction. Each box also contains the sequence 5′-TTTAGTTGGGGTTCAGTTGGCGCGCGCGCGCGC-3′ hybridized to either both (off-state) or the latter sequence (on-state) and one AMP molecule for the on-state. The AMP compound (Figure S2 in Supporting Information)w a s generated with the help of the Maestro Schrodinger program. 35 This software was also used to graft the concerned DNA sequences on the gold surface in combination with the xleap module of AmberTools.

AMP docking computations were performed with the Autodock 4.2 software 36 specifically modified with the help of an artificial intelligence procedure to cope with nucleic-acid targets. 36 This software modification has already been validated and proven its reliability. [37][38][39] The Lamarckian Genetic Algorithm was employed 40 for the global optimum binding position search. For each DNA aptamer conformers, grid boxes enough to encompass the whole target structures were created around the DNA conformers so that the AMP ligand may interact in any target location without a priori. A total of 100 docking poses were realized for each target, and the resulting ligand structures were then clustered according to a root-meansquare deviation (rmsd) lower than 2 Å. The conformation selected was the one which presented the best docking score in the most populated cluster. 40 2.3. MD Simulations. All MD simulations were realized with the Amber 20 30 software suite under periodic boundary conditions with the particle mesh Ewald method. 41 The systems were solvated using TIP3P water molecules 42 with a total of 43,012 molecules for the off-state and 42,749 molecules for the on-state. Ionic concentrations were adjusted to 160 mM with NaCl in order to fit usual experimental conditions, this value being based on the Na + concentration while Cl-ions were added up to the system's neutrality in electric charge, leading to the addition of 129 Na + ions and 60 Cl -ions for both systems. The OL15 force field 43 was used for DNA, and AMP was considered with the general amber force field (gaff2). 44 The OL15 force field was chosen due to the positive results found by Galindo-Murillo et al. for the study of double-stranded DNA when coupled with TIP3P water molecules. 45 The atomic partial charges of the AMP molecule were derived through the AM1-BCC method 46,47 as its reliability has been demonstrated. 48 Consideration of the gold surface functionalized with a 6-mercapto-1-hexanol selfassembled monolayer were reproduced from previous works of the team. 24,33 MD simulations started with a first minimization before engaging a MD heating phase of 150 ps in order to reach 300 K; the last temperature was chosen for all MD simulations. A MD in NTP conditions was then engaged for 500 ps in order to equilibrate the solvent density. The final zdimension lengths of the cells after this process were 16.85 nm for the off-state cell and 15.44 nm for the on-state cell. The simulation systems were then switched to the NVT ensemble, more suitable for simulations with a solid/liquid interface. 15 MD productions (50-200 ns) were made with the help of the SHAKE algorithm 49 in order to use a time increment of 2 fs.

Other simulations were performed on both on-and off-state systems as well as on the aptamer hairpin in isotropic conditions with a similar protocol but staying in the NTP ensemble during the production step of 100 ns.

2.4. Trajectory Analyses. VMD software 50 was used for trajectory visualization and figure generation while other analyses were performed with the cpptraj module of AmberTools. 51 Free energies of binding between the AMP ligand and its aptamer were evaluated with the MMGBSA technique, 52 for which we already demonstrated its reliability in various contexts of ligand associations. [53][54][55] These last computations were made over the last 50 ns of trajectories. Free energy landscape was deduced from trajectory analyses by using Boltzmann statistical mechanics, as previously reported by Ruan and collaborators. 24 3. RESULTS AND DISCUSSION 3.1. Construction of the Off-State System. Inspired by the work made by Yu and collaborators, 28 we designed the 70base DNA three-strand system presented on Scheme 1a. The structure is composed of a central 33 mer sequence containing the anti-AMP sequence, hybridized by 13 base pairs on each extremity to two other sequences 18-and 19-base long, which form the two legs grafted to the surface through a C 6 H 12 -SH function. One key element of our study is the analysis of the capacity of this nucleic acid to achieve the expected arch organization after we graft the thiolated extremities to the surface. We desire that the nucleic acid structure attains a stable U-shaped folding, where both extremities point toward the same direction in space, which would translate in a direct grafting on the surface. Consequently, in the following, the U-shaped geometry of the DNA will be a noteworthy element of our structural enquiries since it allows direct immobilization of our nucleic acid onto the gold/MCH surface. Furthermore, we need to be vigilant when defining the grafting points of our off-state system as they can dramatically influence the overall arch structure and therefore our final results.

For this reason, prior to immobilizing this system on a functionalized gold surface, we thus investigated the organization of this nucleic acid structure in isotropic conditions. The triple-strand structure contains distinct areas of various flexible behavior, being composed of two double-helix sections connected to a central single strand of nucleic acid. When we consider that a single strand DNA is overall more flexible than a double helix, it is reasonable to consider that the central part, which bridges the two double-stranded structures, has the most important role on the curvature of the system. Therefore, we first decipher the conformational space of this central section, that is, the seven unpaired bases with four base pairs on each side, forming a 23-base system (see Scheme 1b).

We ran six replicas of MD simulations in the explicit solvent for the 1b system. To elucidate the representative structure of the 1b system and thus the structural organization of the "arch" with its "two legs", that is, the dehybridized versus the doublestranded parts of the nucleic acids, a free energy landscape was computed from the MD trajectory. Two distances, D1 and D2, corresponding to the distance between the phosphorus atoms of the nucleotides C5/C12 and T3/G13, respectively, were employed as coordinates for the free energy map. These two distances are indicated in Figure 2a. Only the last 100 ns of the six MD replicas were used, leading to a statistical view of 600 ns. This resulted in the free energy landscape shown in Figure 2b.

Figure 2 shows a large zone in the free energy landscape that can be easily accessible, showing that the b system is highly flexible. We were mostly interested in the energy minima as they would be representative of the most stable structures. Among the minima, we actively looked for the one that would give us a U-shaped global structure when transposed back to the initial DNA sequence. These data need to be taken only qualitatively due to the fact that this free energy technique is an approximation and depends heavily on the thoroughness of the conformational sampling.

To get back to the structure of the full 70-base system a defined in Scheme 1a, the values for D1 and D2 were then used as restraints (5 kcal•mol -1 •Å -1 ), taken from the minima highlighted in circles in Figure 2. After 20 ns of MD simulations with restraints followed by 50 ns of relaxation, an appropriate U-shaped structure was extracted from the simulation using D1 between 13.50 and 15.50 Å and D2 between 12.00 and 15.00 Å (area circled in red in Figure 2). Illustrations of all structures obtained in this process are available in Supporting Information. It can be pointed out that despite the weak restraints force constant, the whole structure readily adapts to the restraints, emphasizing the stability of this U-shaped structure. This last geometry was then grafted through a phosphate linker to the Au/MCH surface on both extremities, forming the off-state bridge system shown in the left part of Figure 3. In order to build the onstate system bound to AMP, we first elucidated the structure of the aptamer-AMP complex. To date, there is no experimental 3D structure of this aptamer-AMP complex. The only available description corresponds to three sequences of AMP aptamers grafted on magnetic beads and used in a fluorescence-colorimetric strategy. 29 The prediction made with the mfold software of the two-dimensional organization of these three DNA sequences 31 unveil, in all cases, a hairpin organization in accordance with the experimental results found by Song et al. 29 From these results, we selected the sequence shown in Scheme 1 (bases in bold), which according to their The Journal of Physical Chemistry B based on dielectric property change, we analyzed the ionic distribution over the last 100 ns of each replica by counting the ions over slabs of width 4 Å, parallel to the surface. 24 To obtain an average atom count, the total count was divided by the number of analyzed frames.

When a pure MCH monolayer is considered, Na + and Cl - concentration profiles have the same feature and are both homogeneous above a distance to the surface of approximately 2.5 nm (Figure S5 in Supporting Information). In the case of a mixed MCH/DNA monolayer, the Cl -concentration profile shows almost the same feature as for the pure MCH case, whereas the Na + concentration strongly increases close to the surface for both the on-and off-state (see Figure 6a).

When comparing Na + concentration profiles between both states, a noticeable difference is observed in the ion concentration profile. The peaks on the concentration profile of Na+ extend up to ca. 10 and 8 nm from the surface for the on-and the off-state, respectively. There is a noticeable correlation between the localization profiles of Na + and phosphates in both states (Figure 6b), indicating that the ionic distribution change is directly connected to the architecture conformation switch upon the AMP binding.

CONCLUSIONS

We have proposed a novel DNA-based aptasensor structure for AMP grafted on an Au/MCH surface and we have studied a proposed structure by molecular modeling. We have performed a conformational study on the DNA system in solution in order to optimize its structure before grafting it to the surface to build our initial structure without AMP (offstate). We have elucidated the structure of an aptamer-AMP complex and used it to build our system while in the presence of AMP (on-state). The expected conformational switch in the presence of AMP presents a change in ionic distribution at the vicinity of the surface (around 10 nm) that might be detected by methods based on dielectric property change such as electrochemical and field-effect transistors (capacitance change) and also optical devices (change in the refractive index). This computational study opens the possibility to modulate, in silico, the sequence of each partners in order to screen many different architectures and drive aptasensor design.

As a prospect, and with the will to be as close as possible of real biosensors, especially transistors and electrochemical sensors, we plan to study the effects of electric charges present over the surface on the whole DNA structure and by consequence on the ionic distribution close to the surface.

■ ASSOCIATED CONTENT

* sı Supporting Information

The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jpcb.1c05341. Details on the molecular docking results, schematics for the transferring of the parameters found in the 1b system to the 1a system and the computational pipeline for the prediction of the aptamer/AMP complex structure, snapshots of the structures found after applying the restraints based on the energy minima found after applying the restraints indicated in Figure 2, molecular structure of AMP, rmsd values for one of the on-and offstate simulations each, and ionic profile of a NaCl 160

MM -

  Molecular Mechanics MD -Molecular Dynamics EGOFET -Electrolyte-gated Organic Field Effect Transistor MCH -6-mercapto-1-hexanol AMP -Ampicillin SAM -Self-assembled monolayer VEGF -Vascular endothelial growth factor REMD -Replica Exchange Molecular Dynamics GAMD -Gaussian Accelerated Molecular Dynamics RMSD -Root-mean-square deviation RMSF -Root-mean-square fluctuation LJ -Lennard-Jones SASA -Solvent-Accessible Surface Area PBC -Periodic Boundary Conditions OSC -Organic Semiconductor US -Umbrella Sampling PMF -Potential of Mean Force Chapter 1: Introduction and Bibliographic Review 1.

Figure 1 .

 1 Figure 1. Generic structure of a biosensor.

Figure 2 .

 2 Figure 2. Reproduced from Guo et al 46 . General functioning of their COVID-19 sensor.

Figure 3 .

 3 Figure 3. General structure of an EGOFET biosensor.

  by analyzing the interface through the Gouy-Chapman-Stern model for the electric double layer, by arguing that even though the interactions take place beyond the Debye's length, they still influence the diffuse layer, therefore also influencing the total capacitance of the interface. In the Gouy-Chapman-Stern model, there are two layers at the vicinity of the interface where there is a compact layer composed by solvent molecules followed by a larger diffuse layer in which the ionic concentration decays exponentially until reaching the bulk concentration. An illustration of this model is shown on Figure4.

Figure 4 .

 4 Figure 4. The Gouy-Chapman-Stern model for the electric double layer. In blue, the region containing the rigid Helmholtz layers, in green, the diffuse layer region. In red, the theoretical ionic distribution profile along the distance to the electrode and in purple, the electric surface potential.

Figure 5 .

 5 Figure 5. A typical output curve with five 𝑉 𝐺 values, showing the linear regime in the vicinity of 𝑉 𝐷 = 0.0 V and saturated regime in higher 𝑉 𝐷 values.

  by using a mixed layer of anti-immunoglobulin-G and bovine serum albumin, as seen on Figure6, extracted from their publication. The authors theorize that the transduction process is done by a change in the dipole moment of the anti-immunoglobulin-G molecule attached to the surface generating a change in the surface capacitance. This sensor also proved to be selective between immunoglobulin-G and immunoglobulin-M.

Figure 6 .

 6 Figure 6. Reproduced from Macchia et al 44 , a) general structure of their EGOFET-based biosensor, b) horizontal cross-section of the biosensor showcasing the source (S), drain (D) and gate (G) disposition, the channel and the electrolyte, c) general structure of the functionalized gate interface containing the bioreceptor anti-immunoglobulin-G (anti-IgG) and bovine serum albumin (BSA), d) structural representation of human immunoglobulin-G and -M (IgG and IgM).

  , using instead a peptide aptamer for the detection of the pro-inflammatory cytokine tumor necrosis factor alpha. To conclude, EGOFETs are recognized devices for biosensing, showing high sensibility and versatility in use. Although the concept of a change in interfacial capacitance generating a change in the output current signal is now well-documented, the mechanisms behind this change in capacitance are not always clear. Although there are reports of EGOFET-type biosensors in literature, including using nucleic acids as bioreceptor, reports on the usage of aptamer sequences to recognize low-molecular weight analytes are still scarce.

Figure 7 .

 7 Figure 7. Left: Generic structure of the nucleotide containing the phosphate group, the pentose sugar (ribose for RNA and deoxyribose for DNA) and the nucleobase. Right: Structure of the canonical nucleobases.

Figure 8 .

 8 Figure 8. Reproduced from Wikipedia 71 . The H-bond interactions involved in Watson-Crick base pairing.

Figure 9 .

 9 Figure 9. Workflow of the SELEX method for aptamer selection.

Figure 10 .

 10 Figure 10. Adapted from Pallu et al 67 , a) and b) Output curves of the reported EGOFET with a direct gating by the DNA hydrogel and adding a drop of PBS respectively. c) Evolution over time of the output drain current after the addition of the target with an aptamer functionalized hydrogel (red) and non-functionalized hydrogel (black).

Figure 11 .

 11 Figure 11. Reproduced from Song et al 91 , the secondary structures of the three anti-AMP sequences found in their work with the common regions marked in dashed lines.

Figure 12 .

 12 Figure 12. Reproduced from Taghdisi et al[START_REF] Taghdisi | An Electrochemical Sensing Platform Based on Ladder-Shaped DNA Structure and Label-Free Aptamer for Ultrasensitive Detection of Ampicillin[END_REF] , the mechanism behind their aptasensor for AMP using a "ladder-shaped" DNA structure.

Figure 13 .

 13 Figure 13. Reproduced from Yu and Lai 93 , the mechanism behind their aptasensor for AMP showcasing a ligand-induced conformational switch.

Figure 14 .

 14 Figure 14. Structural formula of Ampicillin.

Figure 15 .

 15 Figure 15. Reproduced from Thyparambil et al 123 . Representative states for unbound and bound native, N-terminal and C-terminal modified NF04 peptide with ochratoxin A (OTA) and ochratoxin B (OTB), equilibrium probabilities of these states (π), and free energies of binding.

Figure 16 .

 16 Figure 16. Reproduced from Phanchai et al 121 , their proposed mechanism for a colorimetric aptasensor for ochratoxin A (OTA).

Figure 17 .

 17 Figure 17. Adapted from Drew and Walsh 128 , exemplary final configurations for their adsorbed hairpin at 300 K, 350 K and 400 K, showing an increase in the system disorder with the temperature.

Figure 18 .

 18 Figure 18. Reproduced from Badhe et al 129 , their protocol for the in silico development of a SARS-CoV-2 biosensor.

Figure 19 .

 19 Figure 19. Reproduced from Ruan et al 130 , their proposed RNA aptasensor showcasing the aptamerflavin complex grafted to a gold surface treated with a MCH monolayer.

6. 2

 2 Charges on surface: the electrode simulations ............................................................................... 6.3 Solvent equilibration ....................................................................................................................... 6.4 Analyses of ion distributions ........................................................................................................... 7. Bibliographic References .......................................................................................................................

  computational techniques and simulations to describe and explore the chemical properties of systems of interest (solids, molecules, solutions, interfaces, among others) and how they evolve in a certain time frame. At the molecular level, it is traditionally divided in two subdomains for energy computation: Quantum Mechanics (QM) and Molecular Mechanics (MM) while Molecular Dynamics (MD) is the study of the evolution of the systems provided by QM and MM over time.

Figure 1 .

 1 Figure 1. The main contributions of a force field potential energy function.

6 )

 6 In this equation, for every pair of atoms i and j in a set of N atoms, we have what we call a 12-6 Lennard-Jones (LJ) potential for the van der Waals interactions coupled with a coulombic potential for the electrostatic interactions. The parameter 𝜀 𝑖𝑗 represents the depth of the LJ potential well where the energy minimum is located.The parameter 𝑅 𝑖𝑗 0 represents the most favorable distance between the two atoms. The distance between the two atoms at a given moment is represented by the term 𝑟 𝑖𝑗 . The factor to the 12 th power refers to the repulsive van der Waals forces when the atoms are too close while the factor to the 6 th power refers to the attractive van der Waals forces occurring when the atoms are sufficiently close to each other. It is important to notice that the Amber force fields do not take into account the van de Waals interactions of vicinal atoms in the same molecule as well as two atoms bonded to the same atom, since they are took into account in the bonding and angle terms. That means that the intramolecular van der Waals interactions are only calculated when, following the chemical bonds, two atoms are at least three bonds (and two atoms) away from each other, in what is called 1,4 van der Waals interactions.

Figure 2 .

 2 Figure 2. Gold/MCH surface used in simulations.

Figure 3 .

 3 Figure 3. Placement of probe on grid covering the entirety of the receptor during molecular docking.

Figure 4 .

 4 Figure 4. Schematic view of the molecular modelling protocol realized during the Umbrella Sampling protocol.

Figure 5 .

 5 Figure 5. Generic structures of the 3-to 6-site models of the water molecule.

Figure 6 .

 6 Figure 6. The boundary effect changing the forces applied on atoms at the surface.

Figure 7 .

 7 Figure 7. Schematic of the Periodic Boundary Conditions.

Figure 8 .

 8 Figure 8. The thermodynamic cycle used in the calculation of the binding free energy of the formation of the complex Aptamer/AMP, with the steps done in gas-phase and in solution.

Figure 9 .

 9 Figure 9. Schematic for the definition of the solvent accessible surface area (orange) and contact surface (green) of a molecule (blue) by the rolling ball algorithm using a spherical probe (red).

Figure 10 .

 10 Figure 10. Structural formula of 6-mercapto-1-hexanol (MCH)

Figure 11 .

 11 Figure 11. Visualization of the PDB structure obtained from the POSCAR file

Figure 12 .

 12 Figure 12. Final surface used on all MD simulations.

Figure 13 .

 13 Figure 13. EGOFET-like circuit showing the hypothetical capacitor formed by the ensemble Gold (Gate electrode)/Electrolyte/OSC.

Figure 14 .

 14 Figure 14. Schematic of the solvent equilibration process of the systems including the MCH/Au surface.

Figure 15 .

 15 Figure 15. The functioning of the counting script for ions. On the left, a schematic with green circles as ions while on the right separated in batches of width 4 Å, a hypothetical histogram describing the distribution of the green circles on each batch.

Chapter 3 : 1 .

 31 Computational studies of a DNA-based biosensor Proposition of DNA-based biorecognition structure for an EGOFET biosensor and expected behavior We propose a biorecognition system composed of three DNA strands grafted on a gold surface treated with a 6-mercapto-1-hexanol (MCH) self-assembled monolayer (SAM) by two anchoring points. The structure is composed of a central 33mer sequence containing the anti-ampicillin (AMP) sequence at one extremity, hybridized by 13 base pairs on each extremity to two other sequences 18-and 19-base long which form the two legs grafted to the surface through a C6H12-SH function.The choice for grafting the DNA strands on Au comes from our efforts in best representing the experimental methods commonly employed for the immobilization of DNA molecules on surfaces[1][2][3] : Au-S bonds are readily formed when gold is treated with alkyl-thiols and MCH is commonly used to create self-assembled monolayers (SAMs) over gold surfaces4 . DNA molecules modified with C6H12-SH terminal functions are commercially available. The proposed biorecognition would take place by a conformational switch undergone in presence of ampicillin (AMP), where the aptamer sequence would de-hybridize and form a hairpin enveloping the AMP molecule, as

Figure 1 .

 1 Figure 1. General scheme of the functioning of the proposed system, passing from the off-state structure (left) to the on-state structure (right) in the presence of the target molecule (AMP) (green), by the formation of the aptameric hairpin (red).

Figure 2 .

 2 Figure 2. Structural formula of Ampicillin (AMP).

1. 1

 1 On the choice of aptamer sequenceThe anti-AMP aptamer sequence used in this work is the following, obtained from literature 2,5 :5′ TTAGTTGGGGTTCAGTTGG 3′Nowadays in literature, there are three described aptamer sequences for the AMP molecule, all originated from the work of Song and collaborators 5 on which they performed magnetic bead-based SELEX (systematic evolution of ligands by exponential enrichment) and utilized them in a dual fluorescence-colorimetric method with Au nanoparticles. Along with the sequence used in this work (named AMP18 by the authors), they have also found the following sequences: 5′-CACGGCATGGTGGGCGTCGTG-3′ (named AMP4 by the authors) and 5′-GCGGGCGGTTGTATAGCGG-3′ (named AMP17 by the authors). All these sequences have in common a GGT(T) section in the loop section with a GC pairing at the beginning of the stem section of their hypothetical secondary structures, which were hypothesized by the authors to play a part on the molecule recognition. These sections are shown on Figure3, extracted from the original article.

Figure 3 .

 3 Figure 3. Reproduced from Song et al 5 , the secondary structures of the three anti-AMP sequences found in their work with the common regions marked in dashed lines.

Figure 4 .

 4 Figure 4. a): 70-bases DNA three-strands off-state system, (bold) anti-AMP sequence, (italic) bases involved in hybridization, b) Structure of the on-state, with the hairpin region marked in red.

Figure 5 .

 5 Figure 5. Inspired by the work of Munzar et al 6 . Evolution of the system free energy during the de-hybridization of an aptamer (red) from a complementary strand (black) immobilized to a surface in the presence of the ligand (green).

Figure 6 .

 6 Figure 6. Schematic on the change of the overall ionic distribution over a surface when a grafted DNA sequence undergoes a conformational sequence. Due to the negative charges present in the DNA molecule, a region around the sequence is expected to have a higher concentration of positive charges (pink). This region shape is expected to follow any conformational change of the DNA molecule.

Figure 7 .

 7 Figure 7. The Gouy-Chapman-Stern model for the electric double layer. In blue, the region containing the rigid Helmholtz layers, in green, the diffuse layer region. In red, the theoretical ionic distribution profile along the distance to the electrode and in purple, the electric surface potential.

Figure 8 .

 8 Figure 8. Process to obtain the off-state DNA structure.

Figure 9 .

 9 Figure 9. The proposed methods of following the folding of the DNA structure. Left: distance between the adenines at the extremities of the system. Right: Radius of gyration, the radial distance from the center of mass that would have the same moment of inertia as the entire system.

Figure 10 .

 10 Figure 10. Left: Evolution of the distance between the adenine bases at the extremities of the off-state system during the simulation in implicit solvent. The red line represents this same distance in the original structure built by NAB. Right: RMSD values in Å of the backbone atoms of the off-state system in implicit solvent with the first frame as reference.

  while the red line indicates the distance of 145.82 Å found in the structure built by NAB 10 .

Figure 11 .

 11 Figure 11. Evolution of the distance between the adenine bases at the extremities of the offstate system during the simulation in explicit solvent. The red line represents the distance in the original structure.

Figure 12 .

 12 Figure 12. Representative structure of the conformation with the minimal obtained values for the distance between the extremities, with both poly-A tails pointing to each other.

Figure 13 .

 13 Figure 13. RMSD values in Å of the backbone atoms of the off-state system in explicit solvent with the first frame as reference.

Figure 14 .

 14 Figure 14. RMSF values in Å of the backbone atoms of the off-state system in explicit solvent.

Figure 15 .

 15 Figure 15. Reduced DNA structure (in red) on which we performed our conformational study.

Figure 16 .

 16 Figure 16. Restraint energy potential function used to keep the hybridization of the double strand sections.

Figure 17 .

 17 Figure 17. Histogram of the radius of gyration of the reduced DNA sequence per percentage of frames analyzed.

Figure 18 .

 18 Figure 18. Reasoning behind the choice of two inter-base distances as structural parameters in our conformational study.

Figure 19 .

 19 Figure 19. Top left: Free energy map of reduced system with the local minima chosen for further analysis circled in black. D1 and D2 are interatomic distances (in Å) defined in the text. The dotted line indicates the region where a linear relation between D1 and D2 is observed. A, B, C and D: Final structures obtained for off-state system after applying restraints indicated on Figure 16 followed by relaxation.

Figure 20 .

 20 Figure 20. The complete off-state bridge system grafted to the Au/MCH surface, with Na + (pink) and Cl -(purple) ions in solution and the AMP molecule in green. The aptamer sequence is marked in red. For the sake of clarity, water molecules are not displayed.

3 .Figure 21 .

 321 Figure 21. The on-state structure, showcasing the main DNA structure (left) containing the aptamer hairpin (red) and the AMP ligand (green) and the secondary DNA probe (right).

Figure 22 .

 22 Figure 22. Computational protocol employed for the prediction of the most relevant 3D structure of the AMP aptamer: a) Predicted hairpin structure obtained from mfold; b) Overlap of the five 3D DNA structures after prediction with 3dRNAweb and 50 ns of relaxation; c) Docking results of the five structures with the AMP molecule in green.

Figure 23 .

 23 Figure 23. The interaction cavity (light blue) on the aptamer hairpin found by docking the AMP molecule (Carbon atoms in green, Nitrogen atoms in blue, Oxygen atoms in red, Sulphur atoms in yellow, Hydrogens in grey).

Figure 24 .

 24 Figure 24. Aptamer-AMP structures. Top Left: Structure 2; Bottom Left: Detail on Structure 2 showcasing the amide group in configuration trans. Top Right: Structure 4. Bottom Right: Detail on Structure 4 showcasing the amide group in configuration cis. Nucleic acids are displayed in blue whereas the AMP ligand is colored by atom types (Carbon atoms in green, Nitrogen atoms in blue, Oxygen atoms in red, Sulphur atoms in yellow, Hydrogens in grey), its molecular surface is displayed in transparency.

Figure 25 .

 25 Figure 25. "On-state" DNA structure containing the aptamer hairpin.

Figure 26 .

 26 Figure 26. In this plot, there are no signs of structure stabilization during the simulation time.

Figure 26 .

 26 Figure 26. All-atom RMSD values in Å of the main "on-state" strand containing the aptamer hairpin in solution during relaxation step.

Figure 27 .

 27 Figure 27. a) Snapshot of the on-state DNA structure (blue) containing the aptamer hairpin (red)

Figure 28 .

 28 Figure 28. All-atom RMSD values in Å of the secondary "on-state" probe strand in solution during relaxation step.

Figure 29 .

 29 Figure 29. The complete on-state system grafted to the Au/MCH surface, with Na + (pink) and Cl - (purple) ions in solution and the AMP molecule in green. The aptamer sequence is marked in red.

Figure 30 .

 30 Figure 30. RMSD values (in Ǻ) (all atoms) for the on-state (red) and the off-state (blue) DNA structures over 100 ns of simulation time with the first frame as reference.

Figure 31 .

 31 Figure 31. a) Snapshot of one of the off-state replicas at t = 0 ns. b) Snapshot of the same off-state replica at t = 100 ns. c) Snapshot of one of the on-state replicas at t = 0 ns d) Snapshot of the same on-state replica at t = 100 ns.

Figure 32 .

 32 Figure 32. Graph of the RMSF values of the on-state DNA structure containing the aptamer hairpin (in Å) obtained per residue, in solution (blue) and grafted to the surface (red).

Figure 33 .

 33 Figure 33. Graph of the RMSF values of the off-state (in Å) obtained per base on the off-state system, in solution (blue) and grafted to the surface (red).

Figure 34 .

 34 Figure 34. Top: Snapshots of the off-state (left) and on-state (right). Bottom: a) Graph of the average Na + ion count versus the distance to the surface in nm. b) The phosphate (ROPO3R') - profile. c) The Cl -profile. The continuous lines represent the off-state while the dotted lines represent the on-state.

Figure 35 .

 35 Figure 35. Ionic profile of both Na + and Cl -over a Au/MCH surface without the DNA structure.

Figure 36 .

 36 Figure 36. The definition of coordination numbers in ions, as the number of oxygen (red) and hydrogen (white) atoms directly surrounding the ions up to the cut-off distance values (in blue) defined as the upper boundary the integration of the radial distribution function (𝑔 𝛼𝛽 (𝑟)) of an atom 𝛽 in the coordination shell a central ion 𝛼.

Figure 37 .

 37 Figure 37. Process of generation of mutant sequences by the deletion of terminal T bases close to the aptamer sequence.

Figure 38 .

 38 Figure 38. Na + profile count of the MD simulation of the mutant on-state structures. In black, the profile of the original sequence. In red, the Na + profile of the mutant with one deleted T base. In green, the Na + profile of the mutant with two deleted T bases.

Figure 39 .

 39 Figure 39. Detail on the hairpin structure in a snapshot of the MD simulation of the mutant with one deleted T base.

Figure 40 .

 40 Figure 40. Desired shapes for an effective EGOFET biosensor with their respective hypothetical Na + distribution profile. To enhance the efficiency of capacitance variation when the analyte binds to the DNA, we suggest that the off-state (left) structure should be as close as possible to the electrode whereas the on-state (right) should extend upright to the electrode surface.

Figure 41 .

 41 Figure 41. Illustration made on the off-state system of C1' atoms of the aptameric strand (orange) with those of its complementary sequence (green), used as collective distance for the reaction coordinate during the US simulation.

Figure 42 .

 42 Figure 42. Illustration of the umbrella sampling window at 40 Å reaction coordinate. The two studied strands are now dissociated.

Figure 43 .

 43 Figure 43. Histogram distribution of the population along the reaction coordinate value.

Figure 44 .

 44 Figure 44. Plot of the reaction coordinate value around the time (frame) recorded for the whole trajectories of umbrella sampling simulation.

Figure 46 .

 46 Figure 46. Structures visualization issued from umbrella sampling simulations for the DNA system in water box. A, B, C, D, E and F represent, respectively, the structure with their collective distance value of 8.8, 13.3, 19.3, 25.05, 31.8 and 40.05 Å.

Figure 47 .

 47 Figure 47. Observation of the off-state dissociation within the umbrella sampling simulation. From left to right, the collective distances have values of 8.8, 17.8, 28.8 and 40.05 Å.

Figure 48 .

 48 Figure 48. Illustration of the ion trapping effect under the surface in periodic image box z+1with Na + (pink) and Cl -(purple) ions in solution. The area in red shows a higher concentration of Na + ions than the green section due to electrostatic interactions. For the sake of clarity, water molecules are not displayed.

Figure 49 .

 49 Figure 49. Snapshots of both off-(left) and on-state (right) at t = 100 ns with total surface electric charge equal to 88 e.

Figure 50 .

 50 Figure 50. Average ionic profiles of both off-(continuous) and on-states (dotted) with a surface charged at 50 e. This was calculated over five replicas, with the error bars indicating the standard deviation of the measurement.

Figure 51 :

 51 Figure 51: Average ionic profiles of both off-(continuous) and on-states (dotted) with surfaces with different values of electric charge applied to them (color-coded, indicated in the figure). Atom count = 0 equates to the atom count found at the region between 12 nm and 20

Figure 52 .

 52 Figure 52. 𝛥𝐼𝑜𝑛𝑠 profile for all the electric charges tested.

Figure 1 .

 1 Figure 1. General scheme of the functioning of the proposed system, passing from an off-state structure (left) to an on-state structure (right) in the presence of the target molecule (AMP) (green) by the formation of the aptameric hairpin (red).

Scheme 1 .

 1 Scheme 1. (a) 70-Base DNA three-strand system; (Bold) Anti-AMP sequence, (Italics) Bases involved in hybridization; (b) Smaller DNA structure formed by the central section of the bridge sequence a; (Underlined) Bases involved in the analyses of interatomic distances performed; (c) Structure of the on-state The Journal of Physical Chemistry B

3. 2 .

 2 Construction of the On-State System. 3.2.1. Structure of the Aptamer-AMP Complex.

Figure 2 .

 2 Figure 2. (a) Structural positioning of D1 and D2 with respect to the 1b sequence. (b) Free energy map of 1b sequence with the local minima tested circled in black and the local minimum where we found our U-shaped structure circled in red. D1 and D2 are interatomic distances defined in the text of the underlined bases in the sequence shown in 1b.

Figure 3 .

 3 Figure 3. Complete off-state bridge system (left) and open on-state system (right) grafted to the Au/MCH surface, with Na + (pink) and Cl -(purple) ions in solution and the AMP molecule in green. The aptamer sequence in both systems is marked in red. For the sake of clarity, water molecules are not displayed.

Figure 5 .

 5 Figure 5. Graph of the RMSF values of the off-state (in Å) obtained per base on the off-state system, in solution (blue) and grafted to the surface (red). Each tick on the x-axis refers to a residue position in the sequence shown below the graph.

Figure 6 .

 6 Figure 6. (a) Graph of the average sodium ion count versus the distance to the surface in nanometers. (b) Graph of the average the phosphate (ROPO3R′) -ion count vs the distance to the surface in nanometers. (c) Graph of the average chloride ion count versus the distance to the surface in nanometers. In all graphs, the continuous lines represent the off-state, while the dotted lines represent the on-state.

Table of contents 1 .

 of1 Molecular Mechanics ................................................................................................................................ 1.1 General concepts about force fields ............................................................................................... 1.2 The Amber force fields .................................................................................................................... 1.3 Specific force-field parameters for gold and MCH .......................................................................... 2. Molecular docking ................................................................................................................................. 2.1 Searching algorithm ........................................................................................................................ 2.2 Scoring function .............................................................................................................................. 3. Molecular Dynamics .............................................................................................................................. 3.1 Equations and time steps ................................................................................................................ 3.2 Free energy determination: the Umbrella Sampling technique ..................................................... 4. Solvent considerations during MD simulations .................................................................................... 4.1 Explicit solvation ............................................................................................................................. 4.1.1 Water models for explicit solvation ......................................................................................... 4.1.2 Periodic boundary conditions .................................................................................................. 4.2 Implicit solvation ............................................................................................................................. 4.2.1 Generalized Born (GB) model .................................................................................................. Classical trajectory analyses ............................................................................................................ Specificities of the Au/MCH system ...................................................................................................... 6.1 Preparation of the surface ...............................................

	4.2.2 Poisson-Boltzmann (PB) method .............................................................................................
	4.2.3 Non-electrostatic contributions ...............................................................................................
	5. MD trajectory analyses .........................................................................................................................
	5.1 5.2 Binding free energy calculations of complex formation: MM-PBSA and MM-GBSA ......................
	6.

...............................................................

  𝐶 𝑡𝑜𝑟𝑠𝑖𝑜𝑛 𝑁 𝑡𝑜𝑟𝑠𝑖𝑜𝑛 + 𝐶 𝑠𝑜𝑙 ∑ (𝑆 𝑖 𝑉 𝑖 + 𝑆 𝑗 𝑉 𝑗 )𝑒

	𝑖,𝑗	𝐴 𝑖𝑗 𝑟 𝑖𝑗 12 -	𝐵 𝑖𝑗 𝑟 𝑖𝑗 6 ) +	𝐶 ℎ𝑏𝑜𝑛𝑑 ∑ (cos 𝜃) 4 ( 𝑟 𝑖𝑗 𝐶 𝑖𝑗 12 -𝑖,𝑗	𝐷 𝑖𝑗 𝑟 𝑖𝑗 10 ) + 𝐶 𝑒𝑙𝑒𝑐 ∑	𝑖,𝑗	𝑞 𝑖 𝑞 𝑗 𝜀𝑟 𝑖𝑗	+
			𝑖,𝑗	-𝑟 𝑖𝑗 2 ⁄ 2𝜎 2

  of MD simulations use the Taylor's expansion of the function 𝑟 𝑖 ⃗ ⃗ (𝑡) to define

	𝑟 𝑖 ⃗ ⃗ (𝑡 + 𝛥𝑡), that is, the new position of the atom after a time step of 𝛥𝑡. This Taylor's
	expansion is as follows:
	𝑟 𝑖 ⃗ ⃗ (𝑡 + 𝛥𝑡) = 𝑟 𝑖 ⃗ ⃗ (𝑡) + 𝛥𝑡. 𝑣 𝑖 ⃗⃗ ⃗ +	1 2	𝛥𝑡 2 . 𝑎 𝑖 ⃗⃗⃗ (𝑡) +	1 6	𝛥𝑡 3 . 𝑏 𝑖 ⃗⃗⃗ (𝑡) +	1 24	𝛥𝑡 4 . 𝑐 𝑖 ⃗⃗ (𝑡) + ⋯ (9)

Each atom i of mass 𝑚 𝑖 undergoes a force 𝐹 𝑖 ⃗⃗ generated from the gradient of potential energy functional of the force field (∇ ⃗ ⃗ 𝐸 𝑡𝑜𝑡𝑎𝑙 ), therefore the acceleration of a certain atom 𝑎 𝑖 ⃗⃗⃗ is given by the following:

Table 1 :

 1 Geometry and potential parameters of the TIP3P model.

		TIP3P
	Length (OH), Å	0.9572
	Angle HOH, °	104.52
	LJ 12 parameter x 10 -3 ,	582.0
	kcal.Å 12 .mol -1
	LJ 6 parameter,	595.0
	kcal.Å 6 .mol -1
	Charge (O)	-0.834
	Charge (H)	0.417

Table 1 :

 1 molecular docking results.

	Structure	ΔG in kcal/mol	% cluster
	Structure 1	-7.3	22%
	Structure 2	-8.8	27%
	Structure 3	-7.6	28%
	Structure 4	-8.5	30%
	Structure 5	-7.8	55%

Table 2 :

 2 Obtained coordination numbers for the Na + and Cl -ions.

		PBE-D	Off-state		On-state
		Bulk 27	Diffuse	Bulk	Diffuse	Bulk
	Na-O	5.8	5.85	5.99	5.98	5.96
	Na-H	15.3	14.31	15.14	14.56	15.08
	Cl-O	6.2	7.66	7.68	7.70	7.69
	Cl-H	5.1	7.14	7.21	7.23	7.22

Table 3 :

 3 Correlation between added charges on the surface and capacitances.In experimental situations, this disparity in capacitances is related to the quality of the MCH SAM grafted to the surface, with poorer coverages of the Au surface increasing the capacitance of the interface. Even though we did not simulate any partially grafted surface with uncovered sections, we hope that our results can indirectly show how those would behave due to defects in the SAM structure.

	Total charge on surface Capacitance
	13 e	5 uF.cm -2
	25 e	10 uF.cm -2
	38 e	15 uF.cm -2
	50 e	20 uF.cm -2
	63 e	25 uF.cm -2
	75 e	30 uF.cm -2
	88 e	35 uF.cm -2

Table 4 :

 4 Obtained coordination numbers for the Na + and Cl -ions on a charged surface.

		PBE-D	Off-state	77 e
		Bulk 27	Diffuse	Bulk
	Na-O	5.8	5.96	6.05
	Na-H	15.3	14.39	14.67
	Cl-O	6.2	7.68	7.69
	Cl-H	5.1	7.24	7.25

  Les biocapteurs sont des dispositifs analytiques qui combinent la spécificité d'un composant biologique ou biorécepteur, tel que des protéines, des acides nucléiques, des ADN, dans le but de mieux comprendre le processus de transduction en analysant les changements de distribution ionique dans notre système sans ou avec l'AMP. ces étapes conduit à la sélection d'un oligonucléotide qui a une affinité de liaison élevée avec la molécule cible. À ce jour, la seule utilisation efficace de séquences d'acide nucléique d'aptamère dans les EGOFET pour détecter sa molécule d'intérêt a été rapportée par Pallu et al 2 , en utilisant un hydrogel d'ADN contenant une séquence sensible au facteur de croissance endothélial vasculaire VEGF (un biomarqueur du cancer). Cela a été fait en utilisant l'hydrogel d'ADN directement comme électrolyte solide couvrant le transistor. Au final, les auteurs ont estimé un total de 175 ± 22 pg de la molécule cible liée à environ 50 mg d'hydrogel d'ADN. En ce qui concerne la littérature plus étroitement liée à ce travail, dans le même groupe que cette thèse a été réalisée, Ruan et al 3 ont publié leurs travaux sur un aptacapteur à ARN Cette thèse peut être considérée comme une suite et un approfondissement des travaux de Ruan et al, dans lesquels les méthodes initialement développées sont poussées plus loin afin de proposer un mécanisme de transduction pour un aptacapteur ADN EGOFET. La chimie computationnelle est le domaine de la chimie sur lequel nous utilisons des techniques de calcul et des simulations pour décrire et explorer les propriétés chimiques des systèmes d'intérêt (solides, molécules, solutions, interfaces, entre autres) et comment ils évoluent dans un certain laps de temps. Au niveau moléculaire, elle est traditionnellement divisée en deux sous-domaines pour le calcul de l'énergie : la mécanique quantique (QM) et la mécanique moléculaire (MM) tandis que la dynamique moléculaire (MD) est l'étude de l'évolution des systèmes fournis par QM et MM au fil du temps.Lors de l'étude de systèmes plus grands à des échelles de temps plus longues, l'approche privilégiée pour définir les interactions entre les particules est la mécanique moléculaire (MM). En MM, nous supposons la validité de l'approximation de Born-Oppenheimer, qui dicte que les fonctions d'onde des électrons et des noyaux atomiques peuvent être traitées séparément en raison de l'énorme différence de masse entre les noyaux et les électrons, donc les électrons ne sont pas pris en compte et la charge est définie directement aux centres des atomes. Toutes les interactions interatomiques sont définies par des fonctions d'énergie potentielle appropriées qui sont collectivement appelées « champ de et en présence d'ampicilline (AMP), respectivement. Nous avons effectivement trouvé des conformations stables pour les deux états, y compris l'élucidation de la structure complexe AMP-aptamère par docking moléculaire. Nous avons évalué les effets de la charge de surface dans notre structure d'ADN et nos distributions ioniques : les structures d'ADN sont stables même dans des environnements électrostatiques extrêmes, et à part une forte augmentation de la concentration en Na+ près de la surface, les distributions ioniques sont toujours distinctes entre les états on et off.

	Le gating (application d'une tension sur le terminal d'électrode de grille) des EGOFET
	conduit à la formation de deux doubles couches électriques, une à l'interface entre la grille et
	tissus et même des organites ou des cellules entières, à un composant physique transducteur capable de mesurer et de quantifier l'interaction entre le composant biologique et l'analyte. Pour une utilisation simple d'un biocapteur, y compris par des personnes sans entrainement, il est préférable que la détection se fasse sans ajout de réactifs supplémentaires ni modification de l'échantillon d'origine (détection sans marqueur et sans réactif). Une approche très puissante consiste à utiliser des surfaces modifiées par des biorécepteurs dont les propriétés physicochimiques dépendent de l'environnement voisin. Dans de telles conditions, la formation d'un complexe biorécepteur/cible entraîne un changement direct des propriétés des surfaces. Parmi les techniques capables de détecter des changements sur l'environnement des surfaces modifiées, il y a les transistors à effet de champ organique électrolytique, ou EGOFET, qui ont été efficacement développés comme transducteurs dans la littérature, étant capables de surveiller directement les processus de bioreconnaissance entre un analyte et son biorécepteur à l'une de ses interfaces de son électrolyte. Dans ce travail, nous effectuons des études computationnelles sur un biorécepteur à base d'ADN pour un biocapteur de type EGOFET en utilisant différentes méthodes de calcul. Nous avons choisi l'ampicilline comme molécule d'intérêt pour le biocapteur. L'ampicilline (AMP) est un antibiotique β-lactame qui a été largement utilisé pour traiter les infections bactériennes chez les humains et les animaux depuis les années 1960 grâce à son efficacité contre les bactéries gram-négatives, ce dont les traitements à la pénicilline n'étaient pas greffés sur une surface d'or traitée avec une monocouche auto-assemblée (SAM) de 6-mercapto-1-hexanol (MCH) par deux points d'ancrage. La structure est composée d'une séquence centrale de 33-mères contenant la séquence anti-AMP à une extrémité, hybridée par 13 paires de bases à chaque extrémité à deux autres séquences longues de 18 et 19 bases qui forment les deux pattes greffées à la surface à travers une fonction C6 H12 -SH. Pour un semi-conducteur de type p (exemples : PEDOT-PSS, DPP-DTT), un VG négatif provoque une augmentation de porteurs de charge à proximité de la surface du semi -conducteur tandis que un VG positif épuise la surface des porteurs de charge, alors que pour un système contenant un semi-conducteur de type n (exemples : diimides naphtalène tétracarboxyliques, dérivés de l'azulène) l'inverse est vrai. Cependant, dans les deux cas, le l'électrolyte et une autre à l'interface entre l'électrolyte et le semi-conducteur. Les deux couches ont une capacité élevée (même si l'interface grille/électrolyte est souvent plus élevée); par conséquent, il est possible de faire fonctionner des systèmes EGOFET en utilisant de faibles tensions, avec moins de 1 V. Comme les courants résultants dépendent également de la capacité de ces deux doubles couches, les EGOFET sont intrinsèquement sensibles aux changements de l'environnement électrique dans ces couches. Cela peut être exploité pour créer des biocapteurs à base d'EGOFET via la fonctionnalisation de l'une des doubles couches qui sont assez importants pour impacter leur capacité peut être amplifié grâce à l'effet de champ, c'est-à-dire qu'un petit changement de capacité conduit à une grande variation de courant de sortie. Les aptamères sont des séquences d'acides nucléiques conçues pour atteindre une affinité élevée pour un analyte cible. Les séquences d'aptamères peuvent être obtenues expérimentalement par une méthode appelée SELEX (Systematic Evolution of Ligands by EXponential enrichment). Cette méthode consiste en une synthèse initiale d'un pool aléatoire de séquences d'oligonucléotides jusqu'à 10 15 séquences différentes, suivie d'un processus de sélection dans lequel le pool aléatoire est incubé avec des molécules cibles immobilisées. Dans cette thèse, le substrat de choix est le même et la distribution ionique sur la surface joue un rôle majeur dans le mécanisme proposé pour un changement de capacité sur la surface. force ». Fondamentalement, un champ de force consiste en une série de fonctions d'énergie potentielle qui décrivent toutes les interactions entre tous les atomes d'un système, soit par l'énergie liante (par rapport aux liaisons chimiques) soit par l'énergie non-liante. Les potentiels sont définis à partir de données expérimentales ou de calculs ab-initio localisés et sont appliqués à l'aide de la mécanique classique. La MM est principalement utilisée pour étudier des processus impliquant un plus grand nombre d'atomes ou de molécules, tels que les interactions de surface, la dynamique des biomolécules et les transitions de phase. Dans notre travail, nous avons utilisé les champs de force classiques Amber OL15 pour les structures d'ADN et le champ de force générique GAFF2 pour les molécules organiques (MCH et AMP). Nous montrons également les résultats des études de calcul que nous avons effectuées dans le système : • Le système est suffisamment robuste pour les échelles de temps de mesure et n'a montré aucun signe de dénaturation au cours des périodes de simulation, y compris la déshybridation des doubles hélices d'ADN ni la déformation du du hairpin pendant les simulations. • Nous avons testé les effets des mutations par délétion dans la distribution ionique et nous avons effectué une petite optimisation sur la longueur des séquences. Nous avons également analysé par Umbrella Sampling l'énergie libre de déshybridation de la séquence d'aptamère du reste de la structure du état off et l'avons comparée à la déshybridation du système en solution, en montrant les changements dans le PMF entre le système greffé et le système en solution, ce qui indique une légère modification du mécanisme de déshybridation provoquée par le greffage. Nous concluons que notre système proposé est prometteur pour tout type de transduction électrique basée sur un changement de capacité ou d'effet de champ, et nous avons élucidé certains aspects du processus de bioreconnaissance qui ne seraient pas si faciles à réaliser expérimentalement, tels que les profils ioniques et les mécanismes de déshybridation capables jusque-là Nous proposons un système de bioreconnaissance composé de trois brins d'ADN interfaces du transistor (grille ou canal). Tout changement de la distribution ionique dans les hairpin de l'aptamère. La molécule AMP reste également en place à l'intérieur possibles.

1

. Cette utilisation intensive a fait de l'ampicilline un polluant environnemental important au fil du temps. Afin de détecter la présence d'AMP avec notre biorécepteur, nous avons ajouté une séquence d'aptamère anti-AMP à notre séquence globale d'La bioreconnaissance proposée aurait lieu par un changement de conformation subi en présence d'AMP, où la séquence d'aptamère se déshybriderait et formerait une épingle à cheveux enveloppant la molécule d'AMP. L'état fermé est la structure attendue lorsqu'il n'est pas en présence d'AMP, que l'on appellera désormais « off-state », tandis que l'état ouvert est la structure attendue en présence d'AMP, ou le « on-state ».

Un EGOFET est composé de trois terminaux (source, drain et grille), d'un semiconducteur organique (OSC) et de l'électrolyte. Lorsqu'une tension est appliquée à la borne de drain (VD) le courant résultant (ID) peut être modulée par l'application d'une tension à la borne de grille (VG). Ceci est obtenu par un effet de champ entre la grille et le semi-conducteur organique, à travers l'électrolyte, qui provoque l'accumulation ou l'épuisement des porteurs de charge dans le semi-conducteur, augmentant ou diminuant ainsi sa conductivité.

VG appliqué sur la grille du système a besoin de dépasser en valeur absolue un certain seuil afin d'allumer correctement le transistor; ce biais minimale est appelé tension de seuil (VT). Comme l'ADN est un polyélectrolyte avec des charges négatives dans sa chaîne, notre hypothèse est que lorsqu'il subit un changement de conformation, il modifie également la distribution ionique globale autour de lui en raison des changements dans le champ électrostatique au voisinage de l'ADN. Par conséquent, lorsque l'ADN est greffé sur une électrode et subit un changement de conformation majeur dans la direction normale à la surface de l'électrode, il devrait changer directement la capacité de l'interface électrode/électrolyte. Les séquences liées plus fortement à la cible sont ensuite extraites et amplifiées par une étape de PCR qui permet des mutations aléatoires sur ces séquences, créant ainsi le pool d'oligonucléotides pour une prochaine génération pour l'étape d'incubation. L'itération de plusieurs de (aptamère + biocapteur) pour la flavine en utilisant des simulations de dynamique moléculaire en 2017. Ils se sont intéressés aux effets du greffage de la séquence aptamère à une surface d'or traitée au 6-mercapto-1-hexanol (MCH) par un spacer de 6 thymidines sur la reconnaissance de la molécule cible. Les auteurs ont démontré que le spacer de thymidine ne perturbait pas la reconnaissance de la flavine, cependant le greffage à la surface a favorisé un changement dans la configuration de l'aptamère, même s'il restait environ à 36 Å de la surface dans une position parallèle. Les auteurs ont conclu que ce changement de configuration, qui indiquait une perte d'affinité pour la cible, était dû à une distribution anisotrope des ions sodium en solution, les ions étaient plus concentrés dans la région entre la surface et l'aptamère que dans la région au-delà de l'aptamère. Dans cette thèse, nous montrons nos résultats sur la construction de notre structure de bioreconnaissance modèle à la fois dans les états off et on, représentant le système en l'absence • Nous avons analysé les profils ioniques à l'interface sur les états on et off et nous avons conclu qu'ils sont effectivement distincts ce qui indiquerait que ces états ont une capacité interfaciale distincte, par conséquent, ces états génèrent deux différents signaux de sortie d'un système de EGOFET que ce soit par la présence ou l'absence d'AMP. Nous avons également suivi la solvatation ionique dans nos systèmes, bien que ces résultats n'aient montré aucune différence entre les couches diffuses et le bulk.

•

https://doi.org/10.1021/acs.jpcb.1c05341 J. Phys. Chem. B 2021, 125, 9499-9506

Acknowledgements

Finally, I wish to thank the referees and members of the defense committee for

In order to maintain a healthy transduction when a voltage is applied to the electrode, this effect must not interfere with the ionic distribution close to the surface to the point that the ionic distributions of the off-and on-states are indistinguishable. This problematic was investigated and will be further discussed on the section 5 of this chapter.

Molecular Dynamics of both states grafted on a Au/MCH surface 4.1 Studies of the DNA structures

We have performed MD simulations for both states in 5 replicas for 100 ns for each replica. We monitored the kinetic, potential, and total energies along the trajectory, as well as the densities, pressures or volumes, and temperatures, none of which showed any unusual behavior.

We analyzed the RMSD of the nucleic acids over the 100 ns in both trajectories to observe global molecular conformational changes. These curves, shown on Figure 30, reach a plateau around 10 Ǻ indicating, this way, that the systems are structurally stable. Considering this low value of RMSD and the size of our DNA structures, it is safe to assume that there were no major distortions on the DNA structures when grafted onto the surface. This is very important, as we need the system to be stable on both off-and on-states during measurements, to guarantee the detection of the presence of the molecule of interest. In a real application, we would consider a large population of structures like ours grafted on a surface and the real output value would depend on the average state of all the structures grafted to the surface. The stability of the system over time validates our choice of structure as it indicates that the system is stable enough so that unwanted transitions between states during measurement are unlikely, theoretically increasing the sensibility of the biosensor by reducing the noise on the output signal.

Potential of mean force and opening mechanism

The appropriate exploration of phase space along the reaction coordinate, as proved on the previous section, let us the possibility to extract the potential of mean force (PMF). This was made with the help of the well-known WHAM software 33 which stands for Weighted Histogram Analysis Method. Figure 45 show the PMFs for both off-state systems: in solution and grafted on the MCH gold surface. A bootstrapping analysis was also performed to provide the Y error bars. These ones highlight the good accuracy of the PMF. Several information can be extracted from these two curves: first, the hybridization free energy is of -19.2 ± 0.2 kcal/mol in solution whereas it is of -18.1 ± 0.2 kcal/mol when this DNA system is grafted on the MCH gold surface. As one can expect, when the system is free in solution, the hybridization free energy is more favorable. However, this variation is reasonably low. Besides, the localization of the minimum is roughly at the same DNA separation distance (the reaction coordinates).

All these facts indicate that the grafting on the surface does not induce steric hindrance so that the DNA segment which contains the aptamer sequence is organized similarly We conclude that our proposed system is promising for any kind of electrical transduction based on a change in capacitance or field effect, and we have elucidated certain aspects of the biorecognition process that would not be so easily done experimentally, such as the ionic profiles and the possible de-hybridization mechanisms. We hope our results instigate further studies on this type of system in the future.

We expect our work to enrich the academic discussion on the molecular modelling of aptasensors, nevertheless, we are aware there are many questions we have left unanswered on the topic. Our optimization of the DNA sequence was quite modest. The DNA is most known for its canonical double helix structure, but the molecule goes way beyond that: triplexes, G-quadruplexes, base modifications, branched DNA… These are all possibilities for the design of a (possibly more effective) DNA-based bioreceptor.

We believe our investigation on the structure of the AMP-aptamer complex was thorough, and we have estimated its binding energy by MM-GBSA. However, there are other methods of binding energy analyses that we did not attempt for this system: GAMD, Umbrella Sampling, Steered Molecular Dynamics… These methods can be interesting for further investigation on this type of system.

We have concluded from experimental results that the aptamer sequence is capable of recognizing the AMP molecule while hybridized by induced-fit. However, we did not attempt to simulate this process, which is a natural progression from our work. Our results on the dehybridization of the aptamer by Umbrella Sampling are very promising but still preliminary for the discussion of the entire transition process between states and the effect of the grafting on the hybridization energy.

It is likely that future studies are going to be able to improve our proposed bioreceptor even further, but we hope this thesis will contribute to the development of a novel aptasensor, to the domain of molecular modelling applied to biosensors and to the progress of science in general.

Graphs of the Δion values per charge applied on surface:

Published article by the team: results, has the best binding interaction when bound to a surface. To obtain the 3D structure of the AMP aptamer, a computational pipeline method 56 was undertaken. An illustration of this pipeline is presented in Figure S3. This recent strategy has been specifically created for single-strand DNA aptamers and validated by reproducing numerous experimental structures. 56 The methodology consists in converting DNA into RNA prior to making explorational folding and then transforming back RNA into DNA.

The first step of this protocol is the generation of a 3D structure from its sequence. We selected the five best-scoring aptamer structures from 3dRNAweb and engaged a MD simulation in a water box for 50 ns to relax these structures. At the end of this step, the five structures tended to converge to a consensual organization. To ensure that these five aptamer s t r u c t u r e sa r ea b l et ob i n dA M P ,m o l e c u l a rd o c k i n g computations were engaged. The interacting site for the AMP molecule was defined as the whole hairpin so that the ligand is free to interact with the entire nucleic acid. Molecular docking energies, compiled in Table S1, cannot discriminate the five structures since their variations are below 2 kcal/mol. 40 One common binding cavity is found for the five aptamer structures even though this ligand could interact with any other region of the DNA structure. On the work of Song et al. 29 the authors have compared several DNA sequences of AMP aptamers and speculated that a common sequence of GGTT in the vicinity of a GC base-pair could be the interacting site of the AMP. In our 3D models, we observe that this pattern belongs to the DNA cavity where the AMP interacts, supporting the relevance of our structures. While the interacting cavity on the hairpin is roughly similar for all aptamer structures, their interacting modes evidence the AMP interacting with this cavity in different orientations on each structure.

To definitely state on the molecular structure of the aptamer, supplemental MD simulations for 100 ns in water boxes were realized for the five complexes. From the five structures, only two kept the AMP in place over 100 ns, both shown in Figure 4. We computed the binding energy using the MM-GBSA method and found two favorable energies of ΔG1=-35.6 ± 4.0 kcal•mol -1 and ΔG2=-41.4 ± 6.2 kcal•mol -1 for structures 1 and 2, respectively. Despite the fact that the free energies are, according to the recovery of uncertainty, of similar values, structure 1 was selected since a supplemental hydrogen bond is formed between the aptamer and AMP, and the conformation of the amide bond in AMP (see Figure S2) adopts a favorable trans conformation, which is not the case in structure 2. For all these reasons, structure 1 was selected for further investigation on the on-state system.

3.2.2. Addition of the AMP-Aptamer Complex to the On-State System. The final AMP-aptamer complex was manually linked to DNA leg (see Scheme 1c). A MD simulation of 100 ns was run, and the resulting structure was attached to the gold/MCH system by a phosphate linker. The other DNA leg which is dehybridized from the bridge was investigated by MD simulations in solution for 100 ns. The resulting most representative structure was then attached to the gold/MCH system by a phosphate linker. The surface with these two attachments constitutes the on-state system (Figure 3, right), which was finally investigated by MD simulations in duplicate.

3.3. MD Simulations of the Full Systems. For both states (on-and off-), we monitored the kinetics, potential, and overall energies along the trajectory, as well as the densities, pressures or volumes, and temperatures. rmsd curves were produced for the nucleic acids, along the 100 ns trajectories, to observe global molecular movements. These curves, available on Supporting Information, reach a plateau indicating, this way, that the systems are structurally stable. In other words, no major distortion on the DNA structures grafted onto the surface was observed.

During MD trajectory visualization, no collapse of the DNA structure onto the surface was observed. This indicates the reliability of the self-assembled monolayer, which is usually used to prevent undesirable adsorption of biomolecules onto gold surfaces. In addition, we found that the DNA base pairs remain associated during the whole simulations, emphasizing that the chosen sequences are accurate. The third observation focused on the AMP which remains bound in the on-state system, demonstrating the high affinity of the aptameric sequence. Another importantr e s u l ti s ,f o rt h eo n -s t a t e simulation, the lack of rehybridization of the DNA strands, reforming the off-state. Taken together, these structural information evidence the self-consistency of our approach.

The root-mean-square fluctuation (RMSF) values were computed for all MD trajectories in order to investigate the variation of flexibility between the systems in solution and grafted onto the MCH gold surface. For the on-state, the fluctuations are less informative to compare since, as expected, the two grafted DNA structures make no specific interaction together (data not shown). For the off-state, comparison of RMSF values in solution and onto the surface led to a decrease of these ones, as illustrated on Figure 5. This is notably the case for the extremities which are directly linked onto the surface. Interestingly, a sharp decrease of flexibility is revealed for the central dehybridized DNA section when the system is grafted to the surface. This demonstrates that the impact of the linkage is not limited to the first nucleic acids residues but extended to the whole sequence and stresses the importance of taking this parameter into consideration when selecting an optimal DNA sequence.

3.4. Analysis of the Changes in Ion Distribution. To further investigate if the simulations are able to provide meaningful information on the detection ability of aptasensors