Keywords: Bioluminescence, Oxyluciférine, Chimie Théorique, QM/MM, Simulation Bioluminescence, Oxyluciferin, Computational Chemistry, QM/MM, MD Simulation

P a g e | 4

Résumé La bioluminescence est un phénomène naturel d'émission de lumière visible par certains organismes vivants. Dans la bioluminescence de la luciole, la lumière jaune-verte émise provient de la dé-excitation de l'émetteur de lumière de son état excité à son état fondamental.

L'émetteur de lumière de la luciole, appelé oxyluciférine, est généré à l'état excité après une oxydation en plusieurs étapes, qui implique une molécule organique (la luciférine) et une enzyme (la luciférase). Au-delà de l'attrait visuel de la bioluminescence de la luciole, ce phénomène naturel inspire et motive une recherche scientifique de grande envergure en raison de ses nombreuses applications. Cependant, la taille et la complexité du système bioluminescent de la luciole compliquent sa compréhension. De nos jours, malgré le grand nombre d'études axées sur la compréhension du processus de bioluminescence de la luciole, certains détails fondamentaux tels que les propriétés photochimiques de l'émetteur de lumière ne sont toujours pas analysés. Ce travail de thèse se concentre principalement sur l'étude de la dernière étape de la réaction bioluminescente de la luciole, c'est-à-dire la transition électronique qui produit l'émission de lumière. Cette thèse analyse et discute la structure chimique de la molécule contribuant à l'émission de lumière d'une part et l'environnement d'autre part afin de mieux nous informer sur le changement de couleur de la bioluminescence de la luciole. À cette fin, on simule l'émetteur de lumière de la luciole, l'oxyluciférine, et certains de ses analogues synthétiques dans différents environnements (dans le vide, dans le solvant et en considérant l'environnement protéique) en utilisant des approches théoriques telles que la méthode hybride couplant la mécanique quantique et la mécanique moléculaire (QM/MM), qui permet l'étude et la simulation de l'ensemble du système bioluminescent de la luciole.

This Ph.D. Thesis was carried out within the theoretical chemistry team of the MSME Laboratory of the Gustave Eiffel University (UGE).

I would like to express my warm thanks to the thesis committee: Dr. Aurélie Perrier-Pineau, Dr.

Adèle Laurent, Dr. Jérémie Léonard, Prof. Nicolas Ferré and Dr. Marco Marazzi, for their insightful comments and suggestions. I greatly appreciated the hard questions and the constructive exchanges during the thesis defense.

I would like to express my sincere gratitude to my supervisor Prof. Isabelle Navizet for the continuous support, guidance and confidence that she has always had in me during the three and a half years' Ph.D. Thesis and Master internship. I am deeply thankful for her kindness, infinite patience and availability (even on weekends or holidays) as well as for taking the time out of her busy schedule to answer all my questions. It was a great pleasure and a privilege having her as my Ph.D. supervisor.

Besides my supervisor, I would like to thank Dr. Romain Berraud-Pache and Dr. Cristina Garcia Iriepa for initiating me to the computational approaches and scientific paper writing. I thank them for their invaluable help in analysing the huge amount of results, but also for the great times spent together.

Many thanks go to: M. Bensifia, R. Panzou, R. Manevy, R. Maskri, M. Cheraki, M. Ponce Vergas, J. Sanz Garcia, M. Sahihi, H. Mouhib, P. Kutudila, A. Borrego-Sánchez for their kindness, moral support and interesting discussions (scientific or not) as well as for the great times spent together.

I would like to express my thanks to Dr. Daniel Roca Sanjuán for having welcomed me in his lab two months and for giving me the opportunity to learn and apply different computational methods.

I give a special thanks to my family, who has always encouraged and supported me: my parents, my sisters and brothers, my wife and my in-laws.

Finally, but by no means least, my thanks go to all theoretical chemistry team for the support, help and guidance during these three years. My thanks also go to the Paris-Est University -Ecole doctorale SIE for providing the funding to attend the national/international conferences and for international mobility grant. I am deeply thankful to everyone who have helped me achieve this Ph.D. Thesis.

Introduction

Raphael Dubois, 1 French pharmacologist, was the first who carried out early work on bioluminescence. He discovered in 1887 that bioluminescence involves three compounds: an oxidized organic compound (luciferin), a specific enzyme (luciferase) and oxygen. Further research reveals that there are a large number of oxidized organic molecules in nature emitting light and each molecule has its specific enzyme. In general, both luciferin and luciferase terms are still used to designate these compounds. By etymology, the word "bioluminescence" comes originally from the Greek "bios" which means "living" and the Latin "lumen" which means "light". Therefore, it is known as a living light. Over 700 bioluminescent species have already been identified and the most fascinating examples of this phenomenon can be mainly found in depths of oceans (e.g., anglerfish) and in some seas which shine with the light of bioluminescent bacteria. Moreover, various terrestrial bioluminescence species can be found in the nature such as bugs (e.g., fireflies) and some mushrooms. Each of these living organisms produce light as the result of chemical reactions. However, the reaction mechanism of bioluminescence system leading to the light emission differs from one species to another, and the light emission colour varies, from blue light to red light, as well. In general, bioluminescence is produced by living organisms for hunting prey, finding mates, defending against predators, etc. Marine species like squids use bacterial bioluminescence for counter-illumination camouflage to protect themselves by adapting to overhead environmental light. Terrestrial species such as firefly uses periodic flashing in its abdomen to attract its mates. Throughout history, bioluminescent organisms have been a target for several areas of research. So far, biologists, chemists and engineers are studying the chemical systems and circumstances involved in bioluminescence to exploit and develop applications that make life easier and safer. For instance, bioluminescence imaging has emerged as a powerful technique for viral infection study. It captures the light emitted from different luciferases to detect viral infection sites and quantify viral replication in small living animal.

Of the known bioluminescent mechanisms, firefly's bioluminescence process is particularly well studied and consists of multistep oxidation. [2][3][4] Firefly's bioluminescence is caused by the oxidation of an organic molecule (luciferin). The D-Luciferin reacts with ATP (adenosine triphosphate) in the presence of a cofactor, Mg 2+ (magnesium ions), inside the cavity of luciferase (biocatalyst enzyme) to form a D-luciferyl-adenylate-luciferase complex. The P a g e | 12 obtained complex is subsequently oxidized by dioxygen and transformed into firefly dioxetanone (also called a high in energy intermediate) and adenosine monophosphate (AMP).

The decarboxylation of firefly dioxetanone yields the oxyluciferin in its first singlet excited state (ES), which decays to the ground state (GS) by emitting a luminous photon in the visible spectral range. The study of bioluminescence system is a great challenge in the experimental field or even in the computational modelling field. Despite the multiple studies focused on understanding the formation and the transformation of firefly dioxetanone (the step 2 and 3 in Figure 1) yielding to the oxyluciferin product as well as the light emission process (last step in Figure 1), several details in bioluminescence system remain unsolved. In this Ph.D. Thesis, I will mainly focus on the last step of the firefly bioluminescent reaction, that is, the electronic transition that yields the light emission (step 4 in Figure 1).

Firefly's bioluminescent colour can be observed in a large ultraviolet-visible spectral window, from yellow-green to red depending on diverse conditions like pH, 5 temperature, [START_REF] Mochizuki | Robust red-emission spectra and yields in firefly bioluminescence against temperature changes[END_REF] luciferin structure modifications [START_REF] Mofford | Aminoluciferins Extend Firefly Luciferase Bioluminescence into the Near-Infrared and Can Be Preferred Substrates over D-Luciferin[END_REF][START_REF] Reddy | Robust Light Emission from Cyclic Alkylaminoluciferin Substrates for Firefly Luciferase[END_REF] and luciferase mutations. [START_REF] Wang | Impact of site-directed mutant luciferase on quantitative Green and Orange/Red Emission Intensities in Firefly Bioluminescence[END_REF] Depending upon conditions, the chemical structure of the light emitter can take six different forms due to inter-exchange reactions such as keto/enol tautomerization and the phenol and/or enol deprotonation/protonation equilibria (Figure 2). This triple equilibrium makes the oxyluciferin forms sensitive to the solvent, pH changes, specific interaction with bases and polarity of the luciferase active site. All of these factors complicate the understanding of the bioluminescence process in particular, the oxyluciferin chemical form contributing to the light emission. A theoretical investigation study on the oxyluciferin chemical forms in different solvents and considering the protein environment has been published in 2011 by Chen et al. [START_REF] Chen | Systematic theoretical investigation on the light emitter of firefly[END_REF] The results of this study exclude all the neutral oxyluciferin forms as light emitters while, anionic oxyluciferin forms are considered as emissive forms. One of the six possible firefly oxyluciferin chemical forms, named phenolate-keto form, has been commonly considered as the most likely form of the light emitter. [START_REF] Chen | Systematic theoretical investigation on the light emitter of firefly[END_REF][START_REF] Naumov | Spectral-Structural Effects of the Keto-Enol-Enolate and Phenol Phenolate Equilibria of Oxyluciferin[END_REF][START_REF] Branchini | Yellow-green and red firefly bioluminescence from 5,5-dimethyloxyluciferin[END_REF][START_REF] Rebarz | Deciphering the protonation and tautomeric equilibria of firefly oxyluciferin by molecular engineering and multivariate curve resolution[END_REF][START_REF] Naumov | Spectral-Structural Effects of the Keto-Enol-Enolate and Phenol Phenolate Equilibria of Oxyluciferin[END_REF] However, there is currently no experimental evidence to confirm the chemical form responsible to the light emission. Despite the numerous theoretical and experimental studies completed throughout the last decade, [START_REF] Navizet | Color-Tuning Mechanism of Firefly Investigated by Multi Configurational Perturbation Method[END_REF][16][17][START_REF] Ugarova | Protein structure and bioluminescent spectra for firefly bioluminescence[END_REF][START_REF] Hirano | Spectroscopic Studies of the Light-Color Modulation Mechanism of Firefly (Beetle) Bioluminescence[END_REF] the debate on this problematic is still ongoing. Indeed the fact that firefly oxyluciferin can exist in six different forms, the bioluminescence emission spectra can then be observed in a large ultraviolet-visible spectral window, from yellow-green to red. [16][17][START_REF] Ugarova | Protein structure and bioluminescent spectra for firefly bioluminescence[END_REF] In this Ph.D. Thesis, I have studied with my collaborators (group of: Dr. Pascal Didier, Dr.

Jérémie Leónard from Strasbourg University and Prof. Nicolas Ferré from Aix-Marseille University, France) the six different chemical forms of oxyluciferin and some of their synthesised analogues in water solution to demonstrate whether oxyluciferin synthetic analogues can be used to mimic the natural oxyluciferin forms in water solution (Figure 3). The analogue forms have been synthesised to avoid the inter-exchange reactions, that is, the triple equilibria. For comparative reasons, the simulations have been performed in explicit water as the realised experimental study by our collaborators was done in water solvent. I present the simulation of the absorption and emission spectra of the chemical compounds in water solution considering the interactions between the compounds and the surrounding solvent molecules by using the P a g e | 14 classical molecular dynamics (MD) and the hybrid quantum mechanics/molecular mechanics (QM/MM) methods. The theoretical and the experimental results were presented in the Journal of Chemical Theory and Computation (JCTC) in 2018. I will present and discuss them in Chapter 6. [START_REF] Garcia-Iriepa | Simulation and Analysis of the Spectroscopic Properties of Oxyluciferin and Its Analogues in Water[END_REF] ). The principal moieties of oxyluciferin structure is given. The phenol-cycle presents a cyclopropyl substituent in the thiazole moiety to block the keto/enol tautomerization. The phenol-OMe presents a methoxy group in the enol group of the thiazole moiety to block both the keto/enol tautomerization and enol deprotonation. While, the OMe-enol presents a methoxy group in the phenol group of the benzothiazole moiety to block the phenol deprotonation.

Thanks to the synthesis of these analogues it became feasible to study experimentally isolated chemical forms of oxyluciferin that could not otherwise be possible. Moreover, it is possible to investigate computationally the photochemical properties of the less stable oxyluciferin forms by using appropriate methods. In addition, computational approaches offer not only the possibility to investigate the firefly's bioluminescence mechanism, but also provide information about the relative stability of oxyluciferin forms, their pKa values and emission colour tuning.

P a g e | 15

Furthermore, I present in Chapter 7 the solvation effects on the absorption and emission spectra of firefly light emitter and its synthetic analogues (Figure 3) by using different solvation models, explicit and implicit models. In addition, I demonstrate which solvation model can accurately reproduce the experimental results. 16 Moreover, I present the simulation of vibrational spectra in different environments as considering the vibronic information in the calculations could be essential to reproduce the overall spectral shape. [START_REF] Jacquemin | Impact of Vibronic Couplings on Perceived Colors: Two Anthraquinones as a Working Example[END_REF][START_REF] Charaf-Eddin | Vibronic Spectra of Organic Electronic Chromophores[END_REF][START_REF] Santoro | Going beyond the vertical approximation with timedependent density functional theory[END_REF][START_REF] Ferrer | Insights for an accurate comparison of computational data to experimental absorption and emission spectra: beyond the vertical transition approximation[END_REF][START_REF] Pedone | Environmental effects in computational spectroscopy: accuracy and interpretation[END_REF][START_REF] Bloino | Aiming at an accurate prediction of vibrational and electronic spectra for mediumto-large molecules: an overview[END_REF] The results of this study were published in the Journal Physical Chemistry Chemical Physics (PCCP) in 2019.

The photochemical properties of the light emitter have been found to be dependent on pH value. 5 The change of the pH value can also modify the oxyluciferin surrounding, for instance the protonation state of both the protein active site and the adenosine monophosphate (AMP).

The change of the protonation state of the AMP inside the protein cavity has been investigated by using the fragment molecular orbital calculations. 27 According to the results of this investigation, the change of AMP protonation state can be responsible for the phenolate-keto polarisation in protein active site. Recently, the interplay between the protonation states of oxyluciferin-AMP complex inside the protein cavity has been studied by using QM/MM methods and it has been found that the AMP deprotonation leads to a blue-shift of the emission spectrum. [START_REF] Garcia-Iriepa | Effect of Protein Conformation and AMP Protonation State on Fireflies' Bioluminescent Emission[END_REF] AMP can exist in two protonation states, mono-and doubly-deprotonated, and hence it can form different H-bonding with the neutral oxyluciferin analogue and its deprotonated form. Then, the change of its protonation state could affect the chemical form of the oxyluciferin and therefore modify the absorption energies. To demonstrate whether the presence of the AMP at different pH values modifies the absorption spectra of the oxyluciferin analogues in water solution, our collaborators (group of: Dr. Pascal Didier from Strasbourg University and Prof. Nicolas Ferré from Aix-Marseille University, France) have evaluated the effect of the AMP protonation state on the absorption spectra of three oxyluciferin analogues (phenol-cycle, phenol-OMe and OMe-enol forms, Figure 3) in water at different pH values in presence or not of AMP. As including the vibronic couplings could be crucial to reproduce the experimental results, I have computed the vibronic absorption spectra of the oxyluciferin analogues-AMP complexes to include the vibronic contributions of the chromophores to the simulated absorption spectra. The results of this study were presented in ChemPhotoChem in 2019. I will present and discuss them in Chapter 8.

Firefly's luciferase and luciferin are widely used for bioluminescence imaging, [START_REF] Woodroofe | Novel Heterocyclic Analogues of Firefly Luciferin[END_REF][START_REF] Greer | Imaging of light emission from the expression of luciferases P a g e | 203 in living cells and organisms: a review[END_REF][START_REF] Kuchimaru | A luciferin analogue generating nearinfrared bioluminescence achieves highly sensitive deep-tissue imaging[END_REF] for instance for detecting tumours in vivo. The efficiency of this technique depends on the brightness of the signal, that is, the emitted bioluminescent light. However, the emitted light from luciferin-luciferase system is yellow-green (around 550-560 nm) in colour which can be P a g e | 16 readily absorbed by the haemoglobin and hence is not able to penetrate biological tissues. To overcome this obstacle and improve the signal detection, several experiments have been done to design a system emitting light colour ranged between red and near-IR spectral window as it has been shown that wavelengths of this spectral range are strong enough to penetrate deeper living tissues. [32][33][START_REF] Weissleder | A clearer vision for in vivo imaging[END_REF] Among the possibilities to tune the light colour of firefly bioluminescence from yellow to near-IR spectral window we can mention the chemical modification of the light emitter structure [START_REF] Woodroofe | Novel Heterocyclic Analogues of Firefly Luciferin[END_REF][START_REF] Ioka | Synthesis of Firefly Luciferin Analogues and Evaluation of the Luminescent Properties[END_REF][START_REF] Ikeda | An allylated firefly luciferin analogue with luciferase specific response in living cells[END_REF] and the mutation of the luciferase. [START_REF] Mofford | Aminoluciferins Extend Firefly Luciferase Bioluminescence into the Near-Infrared and Can Be Preferred Substrates over D-Luciferin[END_REF][START_REF] Wang | Impact of site-directed mutant luciferase on quantitative Green and Orange/Red Emission Intensities in Firefly Bioluminescence[END_REF][START_REF] Jones | Orthogonal Luciferase-Luciferin Pairs for Bioluminescence Imaging[END_REF][START_REF] Branchini | Red-emitting luciferases for bioluminescence reporter and imaging applications[END_REF][START_REF] Branchini | Chemically Modified Firefly Luciferase Is an Efficient Source of Near-Infrared Light[END_REF] Despite the efforts made to clarify the colour tuning phenomenon, many details are still unsolved. [START_REF] Navizet | Color-Tuning Mechanism of Firefly Investigated by Multi Configurational Perturbation Method[END_REF][START_REF] Hosseinkhani | Molecular enigma of multicolor bioluminescence of firefly luciferase[END_REF][START_REF] Iwano | Development of simple firefly luciferin analogs emitting blue, green, red, and near-infrared biological window light[END_REF] In the case of the chemical modification of the light emitter structure, various synthetic analogues have been designed to obtain a red/near-IR emission colour. However, the chemical structure modifications are not an easy task and, in some cases, the emission wavelength of analogues can be observed in blue and in red colours or even outside of the colour spectral window, depending on the substitution made in the light emitter structure. [START_REF] Woodroofe | Novel Heterocyclic Analogues of Firefly Luciferin[END_REF][START_REF] Ioka | Synthesis of Firefly Luciferin Analogues and Evaluation of the Luminescent Properties[END_REF][START_REF] Ikeda | An allylated firefly luciferin analogue with luciferase specific response in living cells[END_REF][START_REF] Iwano | Development of simple firefly luciferin analogs emitting blue, green, red, and near-infrared biological window light[END_REF][START_REF] Zhang | Pyridone Luciferins and Mutant Luciferases for Bioluminescence Imaging[END_REF][START_REF] Mccutcheon | Expedient synthesis of electronically modified luciferins for bioluminescence imaging[END_REF][START_REF] Jathoul | A dual-color far-red to nearinfrared firefly luciferin analogue designed for multiparametric bioluminescence imaging[END_REF] In the Chapter 9 of this Ph.D. Thesis I present the influence of structural modifications within the firefly light emitter structure, oxyluciferin, on the colour tuning of bioluminescence.

For this aim, I have studied the natural firefly oxyluciferin (in its phenolate-keto chemical form) and three of its synthetic analogues that tune the bioluminescence colour from blue to red, called benzothiophene_Oxy, dihydropyrrolone_Oxy and allylbenzothiazole_Oxy. These analogues have been designed by modifying some atoms in the benzothiazole and thiazolone moieties of the oxyluciferin structure. I present the simulation of the absorption and emission spectra of the chemical compounds in protein considering the interactions between the compounds and the surroundings (solvent molecules, protein residues and AMPH) by using the classical MD and the hybrid QM/MM methods. Then, I compare the predicted results to the ones measured experimentally. [START_REF] Woodroofe | Novel Heterocyclic Analogues of Firefly Luciferin[END_REF][START_REF] Ioka | Synthesis of Firefly Luciferin Analogues and Evaluation of the Luminescent Properties[END_REF][START_REF] Ikeda | An allylated firefly luciferin analogue with luciferase specific response in living cells[END_REF] Finally, I suggest a novel oxyluciferin analogue which could be a good candidate to emit light in the near-IR spectral window. The results of this study were published in the Journal Physical Chemistry Chemical Physics (PCCP) in 2020.

This manuscript is divided into three Parts. Part I presents the luminescence types, photoluminescence and bioluminescence phenomena. Part II presents briefly the computational methods that have been used to carry out this work. Part III presents the results obtained during this Ph.D. Thesis. This Part III is divided into five Chapters: i) Chapter 5 presents in detail the computational procedures used to simulate the firefly's bioluminescent system and ii) Chapter 6, Chapter 7, Chapter 8 and Chapter 9 are dedicated to the results discussion of the topics addressed above.

Luminescence types

Luminescence is the process of spontaneous emission of electromagnetic radiation from an electronically excited system. In luminescence phenomenon, the light emission is not resulting from heat energy contrary to incandescence, which is light emitted by a material as a result of heating (e.g., sun, stars and metals in a flame glow by incandescence). There are several types of luminescence, each named according to the excitation mode (source of energy), as shown in Figure 4. Each type of luminescence has its own mechanism leading to light emission and leads to variety of applications in different fields. In the following, only the photoluminescence and bioluminescence are going to be reviewed as they are the phenomena studied in this Ph.D.

Thesis.

P a g e | 19

Photoluminescence

The photoluminescence phenomenon usually involves two main steps and illustrated in the Jablonski diagram (Figure 5): i) Excitation (Absorption), in which the system is excited from the ground (S0) to the excited (S1) state upon irradiation. Thus, the system is temporarily reached a higher energy state. ii) De-excitation, in which the system is de-excited from the S1 to the S0 releasing the absorbed energy by radiative (Emission) or non-radiative transitions, depending on whether a photon is emitted or not during the transition, respectively.

After the excitation of the system, the energy is dissipated by vibrational relaxation (VR). This is a very fast transition between vibrational levels of the same excited electronic state (zigzag line in Figure 5). However, if the vibrational levels of the excited electronic states (S1 and S2) are near in energy, the electron can switch between them, leading to internal conversion (IC). Moreover, the system can return to the ground electronic state S0 by fluorescence, emitting a photon. Furthermore, the electron can change its spin multiplicity during the relaxation process, going for example from S1 to the triplet excited electronic state T1. This is the so-called intersystem crossing (ISC) process. The system can return from the T1 to the S0 by phosphorescence, emitting a photon. Therefore, the luminescence phenomenon can be subclassified according to the emission lifetime: (a) fluorescence with lifetime less than 10 -8 seconds and (b) phosphorescence with lifetime more than 10 -8 seconds. The various processes leading to luminescence are illustrated in Figure 5.

Quantum yield or quantum efficiency

The quantum yield or luminescence efficiency is the ratio of the number of emitted photons to the total number of photons used to excite the molecule. The quantum yield, 𝜙, gives a measure of the efficiency of fluorescing molecule. If 𝜙 ≃ 1 that means that around 100% of the excited molecule returns to its ground state by emitting a photon. In this case, the molecule is considered as a very good luminescent molecule. If 𝜙 ≃ 0 that means that around 100% of the excited molecule decays to its ground state without emitting a photon.

Bioluminescence

Several bioluminescent species have already been identified and the most fascinating examples of this phenomenon can be mainly found in depths of oceans (e.g., anglerfish) and in some seas which shine with the light of bioluminescent bacteria. Moreover, various terrestrial bioluminescence species can be found in the nature such as bugs (e.g., fireflies) and some mushrooms (Figure 6). Each of these living organisms produce light as the result of chemical reactions.

However, the reaction mechanism of bioluminescence system leading to the light emission differs from one species to another, and the light emission colour varies, from blue light to red light, as well. For instance, the system of Vargula higendorfii (marine species) requires only a P a g e | 21 luciferin molecule, luciferase, and dioxygen to produce blue light, whereas the firefly system also depends on ATP to activate the luciferin. Another marine species, Renilla reniformis, associate luciferase with another enzyme (Green Fluorescent Protein) before reacting with its luciferin (called coelenterazine). The result excited product, coelenteramide, transfers the energy to the green fluorescent protein, releasing green light emission (without the green fluorescent protein the luciferin-luciferase Renilla system in vitro releases a blue light).

In general, bioluminescence is produced by living organisms for hunting prey, finding mates, defending against predators, etc. Marine species like squids use bacterial bioluminescence for counter-illumination camouflage to protect themselves by adapting to overhead environmental light. Some Australian sea snails (Clusterwinks) use their shell like a lampshade, filtering the light emitted from naturally glowing cells to deter predators (it shines to avoid being eaten). Other species use bioluminescence as a way to hunt prey. For instance, anglerfish uses an illuminated filament on its head to lure prey. Furthermore, Terrestrial species such as firefly uses periodic flashing in its abdomen to attract its mates. The most bioluminescent species emit only one colour, while others can emit more than one colour like the railroad worm where its body glows green and its head glows red. This glow worm uses its light to warn predators that it is toxic (predators know that eating this worm will cause illness and possible death).

General mechanism of bioluminescence

The production of light by living organisms is carried out in different ways depending on the species, but it has a common pattern. Bioluminescence is a type of chemiluminescence in which the light energy is produced by a chemical reaction, but it takes place inside a living organism. The chemical reaction of bioluminescence requires two unique chemical compounds: luciferin (substrate) and luciferase (an enzyme). In presence of dioxygen, this reaction involves also cofactors such as Mg 2+ (magnesium ions) and ATP (adenosine triphosphate). Luciferin molecule reacts with dioxygen to produce light and luciferase takes a role as a catalyst to speed up the reaction, which is mediated by the cofactors. Thus, bioluminescence results from the oxidation of luciferin by dioxygen catalysed by luciferase in presence of cofactors, leading to an excited state oxyluciferin, which decays to the ground state by emitting a luminous photon (step 1 in Figure 7). Moreover, the energy resulting from the oxidation can be transferred to a fluorescent substance which operates in the reactional mechanism (step 2 in Figure 7). In this case, the stored energy in the complex is used to excite this substance which will serve as a photon emitter. Furthermore, bioluminescence results from the conversion of chemical energy into luminous photons.

Applications of bioluminescence

Throughout history, bioluminescent organisms have been a target for several areas of research. So far, biologists, chemists and engineers are studying the chemical systems and circumstances involved in bioluminescence to exploit and develop applications that make life easier and safer (Figure 8). The technique of ATP detection is one of the oldest applications inspired from bioluminescence. It is the simplest and fastest assay for cell viability. Bacterial bioluminescent genes are widely used in genetic engineering as contaminant biosensors for detecting various toxicants. Vibrio fischeri is a marine bacterium used for monitoring water toxicity. For instance, the presence of pollutants in water can be detected by measuring the amount of light output from the bacterial culture. The decrease of light output is the signal of a possible presence of a contaminant. Furthermore, luciferase systems are also used in genetic engineering like genetic reporters. Bioluminescence imaging has emerged as a powerful technique for viral infection study. It captures the light emitted from different luciferases to detect viral infection sites and quantify viral replication in small living animal.

Today, scientists are coming back toward to bioluminescence phenomenon as a potential form of green energy. They are trying to replace our traditional lamps by glowing trees and buildings. For example, they are developing a biological lighting system that takes the natural properties of bioluminescent symbiotic bacteria of squids.

Luciferins and luciferase

Luciferins belong to different chemical groups (aldehydes, flavins, benzothiazoles, tetrapyrroles, etc.) and thus do not have unique chemical structure (Figure 9). For example, firefly's luciferin (named D-luciferin) is completely different from that of luminescent marine bacteria. The majority of bioluminescent marine species have a luciferin called coelenterazine.

There are other types of luciferins which are rarely studied such as Latia luciferin found in freshwater snail, Latia neritoides. Some bioluminescent species produce luciferin on their own and others absorb it through other species, either as food or in a symbiotic relationship. For instance, squids store the bioluminescent bacteria in their organs and form a symbiotic relationship with them.

P a g e | 24

Figure 9 Few examples of luciferins.

Luciferase structures are constituted from a combination of long chains of organic compounds, amino acid residues. These amino acids share a general formula (Figure 10), which consists of an alpha carbon atom bonded to an amino group (-NH2) a carboxylic acid group (-COOH), hydrogen atom and side chain R (chemical groups). If the R group (side chain) is a hydrogen atom, the amino acid is called glycine, while if it is a methyl group, then the amino acid is called alanine. Thus, the chemical nature of R group is the element that gives the identity of the amino acid. In proteins, twenty common types of amino acids are found, and can be classified into several groups: 1) amino acid with electrically charged side chain (e.g., Arginine and Histidine…), 2) amino acid with polar uncharged side chain (e.g., Serine and Threonine…) and 3) amino acid with hydrophobic side chain (e.g., Valine and Tyrosine…). Other types of amino acid can have R groups with special properties. For instance, R groups can contain a thiol (-SH) groups like Cysteine which can form covalent bonds with another Cysteine. R groups can also form a link with their own amino group and form a ring structure (e.g., Proline). Via a condensation reaction, two or more amino acids combine together with the elimination of a molecule of water to produce a linear chain of amino acids, called a polypeptide.

The covalent bond attaching two amino acids known as a peptide link. The two ends of the peptide chain with free carboxyl (-COOH) and amino (-NH2) groups are known as carboxyl (C-P a g e | 25 terminal) and amino (N-terminal) terminals, respectively. The polypeptide chain apart from the R group is known as the backbone, and the R groups are called side chains. It should be noted that when the peptide links are formed and water molecules are released, we have to use the term "amino acid residues" when discussing proteins.

Luciferase active site

The folding of amino acid chains into 3D structures creates a small pocket in the luciferase, called active site (Figure 11b). The active site is a region of the protein in which a specific substrate binds to it. In addition, it is the region where the chemical reaction takes place.

The active site derives its properties from the (adjacent) amino acids that created it. Moreover, these adjacent amino acids can have side chains (R groups) that are small or large, acidic or basic, hydrophobic or hydrophilic. Therefore, they provide to the active site a very specific size, structural shape, and chemical behaviour. Hence, active site is very sensitive to changes in the protein environment. Among the factors affecting the active site and enzyme function, mention should be made of the temperature and pH. For example, decreasing or increasing the temperature and the pH values outside of a tolerable range can affects the amino acid residues properties and hence, the chemical bonds in the active site and the substrate binding.

Part II

Theoretical Background

This Part II was written using several books such as: "Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory" of A. SZABO and N.S. OSTLUND, [START_REF] Szabo | Modern Quantum Chemistry : Introduction to Advanced Structure Theory[END_REF] "A Guide to Molecular Mechanics and Quantum Chemical Calculations" of J.H. Warren, [START_REF] Warren | A Guide to Molecular Mechanics and Quantum Chemical Calculations[END_REF] "Computational Chemistry and Molecular Modeling: Principles and Applications" of K.I. Ramachandran, G. Deepa and K. Namboori, [START_REF] Ramachandran | Computational Chemistry and Molecular Modeling: Principles and Applications[END_REF] "Computational Chemistry: A Practical Guide for Applying Techniques to Real-World Problems" of D.C. Young, [START_REF] Young | Computational Chemistry: A Practical Guide for Applying Techniques to Real-World Problems[END_REF] "Mathematical Physics in Theoretical Chemistry" of S.M. Blinder, J.E. House, [START_REF] Blinder | Mathematical Physics in Theoretical Chemistry[END_REF] "e-Design: Computer-Aided Engineering Design" of K.H. Chang [START_REF] Chang | e-Design: Computer-Aided Engineering Design[END_REF] and "Biomolecular Simulations: Methods and Protocols" of L. Monticelli, E. Salonen. [START_REF] Monticelli | Biomolecular Simulations: Methods and Protocols[END_REF]

P a g e | 27

Chapter 2

Computational Methods

In this Chapter, the computational electronic structure methods that have been used in this work will be briefly presented. This Chapter is sectioned as follows: We start by introducing the Schrödinger equation, describing ab initio methods from Hartree-Fock to Post Hartree-Fock methods, followed by a treatment of density functional theory (DFT) and Time-dependent density functional theory (TD-DFT). Some basis sets are also briefly presented.

Ab initio methods

Schrödinger equations

In 1926, Erwin Schrödinger 53 developed a new quantum theory which involves partial differential equations, known as Schrödinger equations. The Schrödinger equations exist in two forms: one consisting of time named as time-dependent Schrödinger equation (TDSE) and the other one, in which time factor is neglected and hence named as time-independent Schrödinger equation (TISE). However, these equations cannot be solved exactly for a many-particles system and hence some approximations are needed to find an approximated solution of these equations.

The state of quantum particle in one dimension is represented by its wavefunction 𝜓(𝑥, 𝑡) which comprises two parts, real and imaginary: 𝜓(𝑥, 𝑡) = 𝜓 𝑅 (𝑥, 𝑡) + 𝑖𝜓 𝐼 (𝑥, 𝑡)

Eq. 1

The probability of finding the particle at position 𝑥 is proportional to the square modulus of the wavefunction:

𝑃𝑟𝑜𝑏𝑎𝑏𝑖𝑙𝑖𝑡𝑦 𝑑𝑒𝑛𝑠𝑖𝑡𝑦 = |𝜓(𝑥, 𝑡)| 2 = 𝜓 𝑅 2 + 𝜓 𝐼 2
Eq. 2

P a g e | 28

Given any initial wavefunction 𝜓(𝑥, 0), we can obtain the wavefunction at later times by solving the time-dependent Schrödinger equation (TDSE),

𝐻 ̂𝜓(𝑥, 𝑡) = 𝑖ħ 𝜕 𝜕𝑡 𝜓(𝑥, 𝑡) = - ħ 2 2𝑚 𝑒 𝜕 2 𝜕𝑥 2 𝜓(𝑥, 𝑡) + 𝑉(𝑥)𝜓(𝑥, 𝑡)
Eq. 3

Eq. 3 is a second-order linear differential equation. The Hamiltonian 𝐻 ̂ of the system is given by the expression -

ħ 2 2𝑚 𝑒 𝜕 2
𝜕𝑥 2 + 𝑉(𝑥), where the first term (-

ħ 2 2𝑚 𝑒 𝜕 2
𝜕𝑥 2 ) represents the kinetic energy of the particle, while the second term (𝑉(𝑥)) represents the potential energy of the particle associated with whatever forces are acting on the particle. 𝑖ħ 𝜕 𝜕𝑡 is the energy operator, ħ is Planck's constant and 𝑚 𝑒 is the mass of the particle.

For time-independent problem, the Hamiltonian is not explicitly time-dependent. If the potential energy in the Hamiltonian is assumed to be time-independent, we can separate the time-dependent Schrödinger equation (Eq. 3) into a time-independent form by using mathematical technique such as separation of variables. We assume that the wavefunction can be expressed as a product of a temporal and spatial function,

𝜓(𝑥, 𝑡) = 𝜙(𝑥) 𝑒 -𝑖𝐸𝑡/ħ
Eq. 4

Inserting Eq. 4 into the time-dependent Schrödinger equation we obtain:

- ħ 2 2𝑚 𝑒 𝜕 2 𝜕𝑥 2 𝜙(𝑥) + 𝑉(𝑥)𝜙(𝑥) = 𝐸𝜙(𝑥)
Eq. 5

The operator on the left side of Eq. 5 expresses the Hamiltonian acting on 𝜙(𝑥), which represents the time-independent Schrödinger equation (TISE). This TISE can be written as

𝐻 ̂𝜙(𝑥) = 𝐸𝜙(𝑥) Eq. 6
where 𝐸 is the total energy and 𝜙(𝑥) is the time-independent wavefunction. The timeindependent Schrödinger equation is the equation that should be used when the Hamiltonian operator is not explicitly time dependent and the system does not change with time.

The Hamiltonian for a system consisting of 𝑀 nuclei and 𝑛 electrons described by

Born-Oppenheimer approximation

The Schrödinger equation can be simplified by assuming that the nuclei do not move, as their motion is very slow compared to the high-speed motion of the electrons. This is known as the Born-Oppenheimer approximation. This approximation states that the motions of the electrons and the nuclei can be decoupled as the mass of nuclei are much heavier than electrons and thus, their velocities are negligible. Hence, one can consider the electrons in a molecule to be moving in a field produced by the fixed nuclei. Within the Born-Oppenheimer approximation, the kinetic energy of the nuclei (term 𝑇 ̂𝑁) is neglected, and the repulsion between the nuclei (term 𝑉 ̂𝑁𝑁 ) is considered constant. This leads to an electronic Schrödinger equation:

𝐻 ̂𝑒𝜓 𝑒 = 𝐸 𝑒 𝜓 𝑒 Eq. 9

with

𝐻 ̂𝑒 = 𝑇 ̂𝑒 + 𝑉 ̂𝑒𝑁 + 𝑉 ̂𝑒𝑒 = - 1 2 ∑ 𝛻 𝑖 2 𝑛 𝑖 -∑ ∑ 𝑍 𝐴 𝑟 𝑖𝐴 𝑀 𝐴 𝑛 𝑖 + ∑ ∑ 1 𝑟 𝑖𝑗 𝑛 𝑗>1 𝑛 𝑖 Eq. 10
In Eq. 10, the terms of Eq. 8 describing the kinetic energy of the nuclei is dropped. As nuclearnuclear Coulomb term is a constant, it has to be included to the electronic energy (𝐸 𝑒 ) to calculate the total energy (𝐸 𝑡𝑜𝑡 ) for the system,

𝐸 𝑡𝑜𝑡 = 𝐸 𝑒 + ∑ ∑ 𝑍 𝐴 𝑍 𝐵 𝑅 𝐴𝐵 𝑀 𝐵>𝐴 𝑀 𝐴=1
Eq. 11

It should be noted that the mass of nuclei does not appear in the electronic Schrödinger equation. This method is considered as the most commonly used for the wavefunctions calculations and it is widely implemented in quantum mechanics computer software.

P a g e | 30

Slater determinant

According to the Pauli Exclusion Principle, wavefunctions of electrons are completely anti-symmetric. This states that under exchange of two particle's coordinates and spins the wavefunction changes its sign. In other words, it assumes that the wavefunction of electron must be anti-symmetric with respect to the interchange of the coordinates of any two electrons,

𝜓(1, 2) = -𝜓(2, 1)
Eq. 12

Slater determinant can describe and satisfy the anti-symmetry principle through an appropriate linear combination of Hartree products 𝜓 𝐻𝑃 , which are the non-interacting electron wavefunctions.

𝜓 𝐻𝑃 (1, 2, … , 𝑁) = 𝜒 1 (1) 𝜒 2 (2) … 𝜒 𝑁 (𝑁) Eq. 13
For instance, consider a two-particle case in which the spin orbitals χ1 and χ2 are occupied. If we put the first electron in χ1 and the second in χ2, we will have:

𝜓 1 2 𝐻𝑃 (1, 2) = 𝜒 1 (1) 𝜒 2 (2)
Eq. 14

If we put the first electron in χ2 and the second in χ2, we will have:

𝜓 2 1 𝐻𝑃 (1, 2) = 𝜒 1 (2) 𝜒 2 (1)
Eq. 15

Finally, we can obtain a wavefunction that meets the requirement of the anti-symmetric principle by taking the appropriate linear combination of these Hartree products, 𝜓 1 2 𝐻𝑃 and 𝜓 2 1 𝐻𝑃 .

Thus, the wavefunction of the N-electron system can be described as an anti-symmetrised product of mono-electronic spin orbitals (χi), so-called Slater determinant:

𝜓(1, 2, … , 𝑁) = 1 √𝑁! | | 𝜒 1 (1) 𝜒 2 (1) … 𝜒 𝑁 (1) 𝜒 1 (2) 𝜒 2 (2) … 𝜒 𝑁 (2) 𝜒 1 (𝑁) 𝜒 2 (𝑁) … 𝜒 𝑁 (𝑁) | | Eq. 16
where, the factor determinant equal, which leads to vanish the Slater determinant (the determinant is zero). This is forbidden by the Pauli Exclusion Principle, which states that no more than one electron can occupy a spin orbital. As an orbital cannot contain more than two electrons, the two electrons must have opposing spins (up spin +1/2 and down spin -1/2). It is convenient to use a shorthand form for a normalised Slater determinant, which only presents the diagonal elements of the determinant and includes the normalisation constant:

𝜓(1, 2, … , 𝑁) = |𝜒 1 (1)𝜒 2 (2) … 𝜒 𝑁 (𝑁)〉 Eq. 17
Furthermore, if we always keep the same order of the electron labels, 1, 2, … , 𝑁, then the Eq. 17 can be written as follows:

𝜓(1, 2, … , 𝑁 ) = |𝜒 1 𝜒 2 … 𝜒 𝑁 〉 Eq. 18

Hartree-Fock (HF) approximation

Among the approximate ways to solve the electronic Schrödinger equation, Ĥ 𝑒 𝜓 𝑒 = 𝐸 𝑒 𝜓 𝑒 , the Hartree-Fock (HF) theory constitutes usually the starting point towards more accurate approximations. The HF theory uses the variational theorem to determine an approximate solution of the fundamental electronic state of a system. It starts from using a simplest antisymmetric wavefunction, a single Slater determinant, to describe the ground state of an Nelectron system:

|𝜓 𝐻𝐹 〉 = |𝜒 1 , 𝜒 2 , … , 𝜒 𝑁 〉 Eq. 19
Thus this approximation of the wavefunction ensures a proper description of the electron which obeys the Pauli Exclusion Principle.

Furthermore, according to the variation principle, the best wavefunction is the one that gives the lowest possible energy. If the function is expressed as 𝜓 𝐻𝐹 and the corresponding energy as 𝐸 𝐻𝐹 , then the Schrödinger equation can be expressed as:

〈𝜓 𝐻𝐹 |Ĥ|𝜓 𝐻𝐹 〉 = 〈𝜓 𝐻𝐹 |𝐸 𝐻𝐹 |𝜓 𝐻𝐹 〉 Eq. 20
given the energy expression

𝐸 𝐻𝐹 = 〈𝜓 𝐻𝐹 |Ĥ|𝜓 𝐻𝐹 〉 〈𝜓 𝐻𝐹 |𝜓 𝐻𝐹 〉 Eq. 21
Having the HF wavefunction, 𝜓 𝐻𝐹 , the energy 𝐸 𝐻𝐹 can be easily calculated. Substituting into the energy expression the Slater determinant for the total molecular wavefunction and P a g e | 32

inserting the explicit form of the Hamiltonian operator yields, after simplification of integrals, the energy in terms of the spatial molecular orbitals 𝜓:

𝐸 = 2 ∑ 𝑛 𝑖=1 𝐻 𝑖𝑖 + ∑ 𝑛 𝑖=1 ∑ (2 𝐽 𝑖𝑗 -𝐾 𝑖𝑗 𝑛 𝑗=1 ) Eq. 22
Where 𝐻 𝑖𝑖 is the electronic energy of single electron 𝑖, moving in the field of the nuclei.

Term 𝐽 𝑖𝑗 is called Coulomb integral and represents the electrostatic repulsion between electrons (i.e. between an electron in 𝜓 𝑖 and one in 𝜓 𝑗 ), and 𝐾 𝑖𝑗 is called an exchange integral, related to the anti-symmetric nature of the Slater determinant. It is noteworthy that (2𝐽 -𝐾) terms are related to the electrostatic repulsion between electrons (considered as true Coulombic repulsion).

By minimising the energy in Eq. 22 with respect to the molecular orbitals 𝜓, one can obtain Hartree-Fock (HF) equation, which determines the optimal molecular orbitals leading to the lowest possible energy. The derived HF equation can be expressed in the following form:

𝑓 ̂(𝑖) 𝜓 = 𝜀 𝜓 Eq. 23

Where 𝑓 ̂(𝑖) is the Fock operator of the form

𝑓 ̂(𝑖) = - 1 2 𝛻 𝑖 2 -∑ 𝑍 𝐴 𝑟 𝑖𝐴 𝑀 𝐴=1
+ 𝑣 𝑒𝑓𝑓 (𝑖) Eq. 24

The 𝑣 𝑒𝑓𝑓 term is the effective potential (called also HF potential) which depends on the spin orbitals of the other electrons. The HF potential (𝑣 𝑒𝑓𝑓 ) is the average potential experienced by an electron 𝑖 due to the presence of all the remaining electrons. Hence, The HF equation (Eq. 23) replaces the electron-electron repulsion by an average field (the central idea in the HF approximation). The HF equation (Eq. 23) is non-linear differential equation in the molecular orbitals 𝜓 and must be solved iteratively. The iterative method for solving HF equation is known as the Self-consistent field (SCF) procedure.

Self-Consistent Field (SCF) procedure and Roothan-Hall equations

Self-consistent field (SCF) approach is an iterative procedure which is used to solve the Hartree-Fock equation. The first step in an SCF procedure is to generate an initial guess of molecular orbitals which will be used for a calculation of the average field, that is, the HF potential (𝑣 𝑒𝑓𝑓 ) seen by each electron. Then, a new better set of molecular orbitals is generated,

P a g e | 33
giving in turn a new electronic field. This procedure is repeated until the set of molecular orbitals and the average field no longer change, that is, until reaching the so-called the convergence criterion (SCF converged orbitals).

For obtaining numerical solution, Roothaan [START_REF] Roothaan | New developments in molecular orbital theory[END_REF] and Hall [START_REF] Hall | The molecular orbital theory of chemical valency. viii. a method of calculating ionization potentials[END_REF] have transformed the HF equation (Eq. 23) into a set of algebraic equation by introducing a set of known basis functions (atomic orbitals, AOs) and expanding the unknown molecular orbiatls (MOs) as linear combination of the atomic orbitals (LCAO). In LCAO approximation, each MO is built from a superposition of AOs belonging to the atoms in the molecule. Thus, LCAO approximation refers to construction of a wavefunction from atomic basis functions. This approximation assumes that the total molecular orbitals wavefunction (𝜓 𝑖 ) can be expressed as linear combination of n basis functions (𝜙 µ ):

𝜓 𝑖 = ∑ 𝑐 µ𝑖 𝑛 µ 𝜙 µ Eq. 25
where 𝑐 µ𝑖 are the MO expansion coefficients. Each MO 𝜓 is expanded in term of n basis functions (𝜙 µ ). Usually, the basis functions are centered at the nuclear positions, they are referred to as atomic orbitals (AOs).

Substituting into the HF equation (Eq. 23) the LCAO approximation (Eq. 25) for the MO's 𝜓 yields the Roothaan-Hall equations, which can be expressed compactly as a set of matrix equations:

𝐹𝐶 = 𝑆𝐶𝜀 Eq. 26
where 𝐹 is the Fock matrix, 𝐶 is the matrix of the MO coefficients (𝑐 𝑢𝑖 ). The matrix 𝑆 gives the overlap between the orbitals. The MO energies are given by the diagonalised matrix 𝜀. The Fock matrix 𝐹 is analogous to the Hamiltonian in the Shrödinger equation and its elements are given as:

𝐹 µ𝑣 = 𝐻 µ𝑣 𝑐𝑜𝑟𝑒 + 𝐽 µ𝑣 -𝐾 µ𝑣 Eq. 27
where 𝐻 𝑐𝑜𝑟𝑒 is the core Hamiltonian of the form: The Fock matrix 𝐹 depends on the coefficients C, which makes the Roothan-Hall equation non-linear and therefore must be solved with an iterative approach such as the SCF procedure. The advantage of the Roothan-Hall equation formalism is that the MOs can be expanded into a linear combination of known one-electron functions as shown in equation (Eq. 25).

𝐻 µ𝑣 𝑐𝑜𝑟𝑒 = ∫ 𝜙 µ (𝑟) [- 1 2 𝛻

Basis sets

A basis set in computational chemistry is a set of mathematical functions (called basis functions), linear combinations of which yield molecular orbitals as:

𝜓 𝑖 = 𝑐 1 𝜙 1 + 𝑐 2 𝜙 2 + ⋯ + 𝑐 𝑛 𝜙 𝑛 Eq. 31
where 𝑛 is the size of the basis set. 𝜙 1 , 𝜙 2 , … , 𝜙 𝑛 and 𝑐 1 , 𝑐 2 , … , 𝑐 𝑛 represent the basis functions and the normalisation constants, respectively. We can also express the equation (Eq. 31) as shown in equation (Eq. 25).

The most popular type of basis sets used are Slater type orbitals (STOs) and Gaussian type orbitals (GTOs). The STOs are obtained from a radial part (describes the radial extend of the orbital) and an angular part (describes the shape of the orbital). STOs have a form depending on spherical coordinates:

𝜙 1 (𝛼, 𝑛, 𝑙, 𝑚; 𝑟, 𝜃, 𝜙) = 𝑁𝑟 𝑛-1 𝑒 -𝛼𝑟 𝑌 𝑙,𝑚 (𝜃, 𝜙) Eq. 32
where 𝛼 is the orbital exponent, 𝑛, 𝑙, and 𝑚 represent quantum numbers: principal, angular momentum, and magnetic, respectively. The 𝑟, 𝜃 and 𝜙 represent spherical coordinates, and 𝑌 𝑙,𝑚 represents the angular momentum part (the function describing the shape of the orbital).

𝑁 is the normalisation constant. In modern quantum chemirstry computations, the STO functions are restricted to atomic and diatomic systems because they are extremely slow when evaluating the two-electron integrals. To overcome this difficulty, some modifications to the wavefunction have been made by introducing Gaussian type functions. Gaussians. However, the valence orbital is now split into three sets of contracted Gaussian functions, 3, 1 and 1 primitive Gaussians, respectively. Orbitals with angular momentums can be added to these basis sets and lead to the so-called polarised basis sets. The polarised basis sets such as 6-311G(d,p) add d-functions to carbon atoms and p-functions to hydrogen atoms.

Correlation consistent basis sets, such as cc-pVDZ (correlation consistent valence double zeta), cc-pVTZ (correlation consistent valence triple zeta) and aug-cc-pVDZ (Augmented versions of cc-pVDZ), etc., are specially designed for post-HF calculations in which electron correlation is better included than at HF level. The "p" stands for "polarisation" indicating that polarisation (correlating) functions are included.

P a g e | 36

Post Hartree-Fock methods

The major limitation in the HF method is that it completely neglects the electron correlation (correlations between electrons with same spin). That is, HF includes only the average effect of electron repulsion while the explicit electron-electron interactions are ignored (each electron sees the other electrons as a mean-field). The correlation energy (𝐸 𝑐𝑜𝑟𝑟 ) is defined as the difference between the HF energy (𝐸 𝐻𝐹 ) and the exact energy of the system (𝐸 𝑒𝑥𝑎𝑐𝑡 ) obtained by solving the Schrödinger equation in the employed basis set,

𝐸 𝑐𝑜𝑟𝑟 = 𝐸 𝑒𝑥𝑎𝑐𝑡 -𝐸 𝐻𝐹 Eq. 34
This correlation energy is a negative quantity because the HF energy is an upper bound to the exact energy.

HF approach is a good starting point for more accurate methods, post Hartree-Fock methods, in which the electron correlation effects is included. In the following, some of the calculation approaches begin with a HF calculation and then correct for correlation effects are going to be reviewed. Among of them: Configuration Interaction (CI), Multi-Configurational Self-Consistent Field (MCSCF), and Complete Active-Space Self-Consistent Field (CASSCF) method as well as its second order perturbation approach (CASPT2).

Configuration Interaction (CI)

In the HF approach, all the electrons are located in the lowest molecular orbitals (MOs).

Therefore, the Slater determinant is built using only the orbitals occupied in the electronic ground state while the virtual (unoccupied) orbitals are not exploited. In the Configuration Interaction (CI) approach, a linear combination of Slater determinants rather than one single determinant in HF is applied to approximate the wavefunction. The linear combination of Slater determinants that describe the many-electron wavefunction can be written as: When Slater determinants are included for all possible excitations, the full configuration interaction (FCI) model is obtained, the solutions of which corresponds to the exact Schrödinger wavefunctions in the used basis set. In other words, in FCI, the electronic wavefunction is expressed as an expansion of all possible configurations constructed by considering electronic excitations from a reference zeroth-order wavefunction. In practice, for a large system, a full CI calculation is computationally infeasible due to the huge number of Slater determinants included in the calculation. One possible method for overcoming this inability is to limit the number of excitations considered in the CI expansion like the truncated models, which include only a specified number of excitations. Most commonly, CI methods are truncated at the CISD (only singly and doubly excitations are included) level since they provide quantitative agreement with the full CI energy (generally represent about 95% of the basis set correlation energy in small molecules at their equilibrium geometries). Including a triply excitations, the percentage of the achieved correlation energy can be increased to 96.73% while, the addition of the quadruple excitation is needed for achieving almost 100% of the correlation energy. The truncated approaches are no longer size-consistent nor size-extensive, meaning that their quality degrades for larger molecules.

𝜓 = 𝑐 0 𝜙 0 + ∑

Multi-Configurational Self-Consistent Field (MCSCF)

In the Multi-Configurational Self-Consistent Field (MCSCF) method, the wavefunction is constructed as a linear combination of Slater determinants. The MCSCF method involves the simultaneous optimisation of both the expansion coefficients and orbitals. are optimised. In CASSCF method, the orbitals are distributed into three subspaces. The first subspace represents the inactive space in which all orbitals are doubly occupied (closed shell) in all Slater determinants. In the second subspace, called virtual space, the orbitals are always unoccupied in all Slater determinants. The third subspace is called active space. In this space, a full CI expansion is considered. In addition, the electron and orbitals included in the active space are those that contribute most to the multi-reference character of the system of interest. The organisation of the molecular orbitals spaces is schematised in Figure 12. virtual (unoccupied), active (valence) and inactive (doubly occupied) orbitals. Six active electrons in six orbitals are included in the active space, leading to a CASSCF (6,6).

In CASSCF calculations, the selection of the molecular orbitals for the active space should be performed by the user. A priori knowledge of the chemical nature of the system under study is needed for the proper choice of the active space. The major problem with CASSCF method is that it includes only the non-dynamic electron correlation for active space orbitals and provides no correlation for the electrons outside of the active space. Thus, the dynamic component of electron correlation is poorly treated. CASSCF calculations are usually indicated as CASSCF(𝑛, 𝑚)

to specify the number of electrons, 𝑛, in 𝑚 orbitals including in the active space. Furthermore, with this method it is possible to optimise either a single state (SS-CASSCF) or state-averaged (SA-CASSCF).

P a g e | 39

Complete Active Space Perturbation Theory

If the active space has been adequately chosen, the CASSCF wavefunction can be used as the starting point for more accurate treatments including multireference perturbation theory as it provides reference configurations for the calculation of the dynamic electron correlation.

Multireference second-order perturbation theories have been developed in order to treat both dynamic and non-dynamic (static) electron correlation. Among these perturbation theories, the Complete Active Space Perturbation Theory to the second order (CASPT2, implemented in MOLCAS and MOLPRO programs) that uses the CASSCF wavefunction as the reference zerothorder function. In CASPT2 methods, the Hamiltonian operator (𝐻 ̂ ) can be expressed by the sum of a reference zeroth-order function and a small perturbation term (𝑉 ̂ ),

𝐻 ̂ = 𝐻 ̂0 + 𝜆𝑉 ̂ Eq. 38
Where 𝐻 ̂0 is the HF Hamiltonian, 𝜆 is an arbitrary real parameter (it should be small in order to obtain accurate result). The perturbation should be small in comparison with the Hamiltonian eigenvalues. The Schrödinger equation for a perturbed system can then be written as

𝐻 ̂𝜓𝜆 = (𝐻 ̂0 + 𝜆𝑉 ̂)𝜓 𝜆 = E 𝜆 𝜓 𝜆 Eq. 39
Furthermore, the wavefunction and the energy can be expanded in terms of a Taylor series:

𝜓 𝜆 = ∑ 𝜆 𝑛 𝜓 𝑛 ∞ 𝑛=0 , E 𝜆 = ∑ 𝜆 𝑛 E 𝑛 ∞ 𝑛=0
Eq. 40

In the CASPT2 calculations, the target excited states are characterised in a qualitative manner by the CASSCF wavefunctions. In the CASPT2 method, these wavefunctions are used as the reference states including the static electron correlation to describe the first-order perturbed wavefunction and second-order perturbation energies to further account for the dynamic correlation. There are two theories based on CASPT2 to excited state calculations:

state-specific and multi-state CASPT2 (MS-CASPT2). In the state-specific CASPT2 theory, the perturbation corrections are separately evaluated for each CASSCF reference wavefunction.

However, the state-specific treatment can be inadequate for studying near-degenerate states, in which the strong mixing occurs between the reference state and one or more secondary CASSCF states. This mixing can properly account for by the multi-states CASPT2 theory, which is an extension of CASPT2 based on the quasi-degenerate perturbation theory in conjunction with the multi-partitioning methods. The MS-CASPT2 uses an effective Hamiltonian matrix in the reference state basis under the presence of the perturbation, and its diagonalization yields the mixing of the perturbed states.

P a g e | 40

The main disadvantages of the CASPT2 method reside in: i) the difficulty in choosing the active space (limited to about 15 electrons in 15 orbitals due to the size of the CI calculations),

ii) computational time needed is significantly larger compared to CASSCF calculations and iii) analytical energy gradients are not feasible for medium-large active spaces.

In the following section, we will give a brief review of the basic density functional theory (DFT) and time-dependent DFT formalisms, which are the most promising approaches for describing the dynamics of interacting many-body systems.

Methods based on electronic density 2.2.1 Density Functional Theory (DFT)

In general, ab initio calculations give good and accurate quantitative results for small systems while they become prohibitively expensive with increasing the system size. The advantage of the correlated methods is that they can treat systematically the electron correlation effects. In addition, when all the approximations are made sufficiently small in magnitude, they can converge to an exact solution. On the other hand, the major drawbacks of correlated calculations are due to the enormous amounts of computer CPU time, disk space, memory and high computation cost. Density Functional Theory (DFT) method is an alternative to multi-configurational methods for determining electronic structures. It has become popular in recent years. DFT is considered as the dominant method for the quantum mechanical simulation of medium-large systems. It permits to achieve a relatively high level of accuracy of medium-large systems at a low calculation cost and computer CPU time. The principle behind DFT is simpler than the one of multi-configurational methods. It assumes that the energy of a molecule can be obtained from the electron density instead of a wavefunction. In addition, the electronic energy is expressed in terms of the electron density, which depends uniquely in three spatial and one spin coordinates. Furthermore, DFT is a variational method based on the Hohenberg-Kohn (HK) theorems. The first HK theorem states that the ground state energy of a system depends only on the electron density, that is, it is a functional of electron density.

Therefore, the total ground state energy 𝐸[𝜌] of an electron system can be written explicitly as The first term in the square brackets is the kinetic energy of the non-interacting reference system, the second term is the external potential, the third one is the Hartree potential, and the last one is exchange-correlation potential. 𝜙 and 𝜀 are the Kohn-Sham orbitals and their energy, respectively. Eq. 44 can be simplified in a more compact form by defining an effective potential 𝑣 𝑒𝑓𝑓 :

a
[- 1 2 𝛻 𝑖 2 + 𝑣 𝑒𝑓𝑓 ] 𝜙 𝑖 (𝑟) = 𝜀𝜙 𝑖 Eq. 45
The KS equation (Eq. 45) looks like Hartree-Fock equation which need to be solved selfconsistently.

Following the KS approach, the 𝐹 [𝜌] term in Eq. 41 becomes:

𝐹 [𝜌] = 𝑇 [𝜌] + 𝑉 𝑒𝑒 [𝜌] = 𝑇 𝑠 [𝜌] + 𝐽 [𝜌] + 𝐸 𝑥𝑐 [𝜌] Eq. 46
where 𝑇 𝑠 represents the kinetic energy of a system of non-interacting electrons, 𝐽 represents the Coulomb repulsion of the electronic density with itself and 𝐸 𝑥𝑐 is the exchange-correlation energy. Inserting Eq. 46 into Eq. 41, the electronic energy can be expressed as:

𝐸[𝜌] = ∫ 𝜌(𝑟) 𝑣(𝑟) 𝑑𝑟 + 𝑇 𝑠 [𝜌] + 𝐽 [𝜌] + 𝐸 𝑥𝑐 [𝜌]
Eq. 47

𝐸 𝑥𝑐 [𝜌] is unknown term and needed to be approximated, which is the focus of the next part.

The kinetic energy of a system of non-interacting electrons 𝑇 𝑠 and the electron density 𝜌 are known exactly from the KS orbitals 𝜙 𝑖 :

𝑇 𝑠 (𝜌) = - 1 2 ∑ 𝑁 𝑖 〈𝜙 𝑖 |𝛻 2 |𝜙 𝑖 〉 Eq. 48 𝜌 = ∑ 𝑁 𝑖 |𝜙 𝑖 | 2
Eq. 49

So far, the exact formulation of the exchange correlation term (𝑥𝑐), which includes exchange, correlation and a small kinetic component, does not exist and hence, different approximations have been proposed to solve the exchange-correlation energy. The simplest approximation is the Local Density Approximation (LDA), which assumes that exchangecorrelation energy of the system is equivalent to the one of a homogeneous electron gas with the density. The expression of the exchange-correlation energy of such a homogeneous electron gas is known as

𝐸 𝑥𝑐 𝐿𝐷𝐴 [𝜌] = ∫ 𝜌 𝜀 𝑥𝑐 𝑢𝑛𝑖𝑓 [𝜌] 𝑑𝑟 Eq. 50
where 𝜀 𝑥𝑐 𝑢𝑛𝑖𝑓 is the term of the exchange-correlation energy per electron in a homogeneous electron gas of density, that is, the exchange-correlation energy density is assumed to be that of a uniform electron gas of the same density. The LDA approximation is, in principle, adequate for systems with a low spatial varying electron density. This approximation can be refined by expressing the different energy terms not only as a function of density but also of its gradient, so-called the Generalized Gradient Approximation (GGA), and its kinetic energy density, socalled meta-GGA approximation (a second derivative of the density is also introduced). In the GGA the gradient of the density (𝛻𝜌) is included in the exchange-correlation energy

𝐸 𝑥𝑐 𝐺𝐺𝐴 [𝜌] = ∫ 𝜌 𝜀 𝑥𝑐 𝑢𝑛𝑖𝑓 [𝜌, 𝛻𝜌] 𝑑𝑟 Eq. 51 P a g e | 43
The most popular GGA functionals are PW91 (functional proposed by Perdew and Wang) 56 , PBE (functional proposed by Perdew, Burke and Ernzerhof) 57 and BLYP (combination of the exchange functional by Becke and the correlation functional of Lee, Yang and Parr) 58,59 . These approximations give reliable results for chemical bonds, but mostly fail for van der Waals forces or dispersion interactions. The hybrid functional is another generation functional in which an "exact" exchange from HF theory is introduced. The most widely used are B3LYP 58,59 and PBE0 60 functionals. The formulation that simplifies the exchange-correlation energy of the hybrid functional (𝐸 𝑥𝑐 ℎ𝑦𝑏 ) can be expressed as: Eq. 53

𝐸 𝑥𝑐 ℎ𝑦𝑏 = 𝑏𝐸 𝑥 𝐻𝐹 + (1 -
Here Δ𝐸 𝑥 𝐵88 and 𝛥𝐸 𝑐 𝑃𝑊91 are GGA corrections used to the LDA exchange and correlation energies respectively. 0.2, 0.72 and 0.81 are three adjustable parameters of B3LYP functional formulation. There is another generation functional that combines a short-range density functional exchange with a long-range HF exchange, called range-separated hybrid functional.

The hybrid functionals and range-separated hybrid functionals have been widely used in quantum chemistry, particularly the B3LYP and cam-B3LYP functionals. The B3LYP functional has been mainly used in this Ph.D. Thesis.

Time-dependent density functional theory (TD-DFT)

TD-DFT is a quantum mechanical framework used in chemistry, biology, and materials science to describe the properties and dynamics of interacting electronic many-body systems in the presence of time-dependent (TD) potentials. TD-DFT extends the basic ideas of ground state DFT to the treatment of excitations or more general TD phenomena. The Runge-Gross TD-DFT theorem (exposed in 1984) 61 proves that, for all observable properties of a many-electron system evolving from in a given initial state, there is a one-to-one correspondence, 𝑣 𝑒𝑥𝑡 (𝑟, 𝑡)↔𝜌(𝑟, 𝑡), between the external potential (𝑣 𝑒𝑥𝑡 (𝑟, 𝑡)) and the electronic one-body density (𝜌(𝑟, 𝑡)). The core of TD-DFT is the Runge-Gross theorem, which is an analogue of the P a g e | 44

Hohenberg-Kohn theorem extended to the TD system. The density 𝜌 (𝑟, 𝑡) is the probability of finding at time 𝑡 one electron at position 𝑟 1 , another at position 𝑟 2 , etc., expressed as:

𝜌(𝑟, 𝑡) = 𝑁 ∫ 𝑑 3 𝑟 2 … 𝑑 3 𝑟 𝑁 |Ψ(𝑟 , 𝑟 2 … 𝑟 𝑁 , 𝑡)| 2
Eq. 54

where the density (𝜌(𝑟, 𝑡)) is normalised at all times to the particle number 𝑁. This density is the basic variable in terms of which TD-DFT is formulated. If we know the TD electronic density of a system, evolving from a given initial state, we can determine the external potential, which produce this density. Therefore, with the external potential, the TD Schrödinger equation can be solved and hence all the properties of the system can be obtained.

Having the Runge-Gross theorem, a time-dependent Kohn-Sham (TD-KS) system can be constructed by introducing an auxiliary system of non-interacting electrons subject to an external potential (𝑣 𝐾𝑆 ). In this theorem, the external potential is chosen such that the density of the KS electrons is the same as the one of the interacting system. The KS electrons obey the following TD Schrödinger equation:

𝑖 𝜕 𝜕𝑡 𝜙 𝑖 (𝑟, 𝑡) = 𝐻 ̂𝐾𝑆 (𝑟, 𝑡)𝜙 𝑖 (𝑟, 𝑡) = [- 1 2 𝛻 𝑖 2 + 𝑣 𝐾𝑆 [𝜌](𝑟, 𝑡)] 𝜙 𝑖 (𝑟, 𝑡) Eq. 55
where 𝐻 ̂𝐾𝑆 is the Kohn-Sham Hamiltonian. The density of the interacting system can then be reproduced from the TD-KS orbitals, that is:

𝜌(𝑟, 𝑡) = ∑ 𝑁 𝑖=1 |𝜙 𝑖 (𝑟, 𝑡)| 2
Eq. 56

where the density 𝜌(𝑟, 𝑡) is the same that one would have determined from the equation (Eq. 54). Therefore, if we knew the exact KS potential (𝑣 𝐾𝑆 ), we would determine from the equation (Eq. 55) the exact KS orbitals and hence, the correct density of the system. The KS system is not intended to reproduce the full many-body wavefunction of the interacting system, only its one-body density. Furthermore, the TD-KS equation features an effective one-body potential, which ensure that the non-interacting system reproduce the correct density. It is usually written in the following form: It uses a simple model of the intra-and inter-molecular interactions within a molecular system. The energy of the system is divided into contributions like the bonds stretching, bond bending, torsions and non-bonding (van der Waals forces and electrostatic interactions) terms.

𝑣 𝐾𝑆 [𝜌](𝑟, 𝑡) = 𝑣 𝑒𝑥𝑡 (𝑟, 𝑡) + ∫ 𝑑 3 𝑟 ′ 𝜌(𝑟 ′ , 𝑡) |𝑟 -𝑟 ′ | + 𝑣 𝑥𝑐 [𝜌](𝑟,
The form of the algebraic expression for the energy, and parameters in it, constitute a force P a g e | 47

field, which contains these four contributions. Thus, the typical form of a force field can be written as: [START_REF] Warren | A Guide to Molecular Mechanics and Quantum Chemical Calculations[END_REF][START_REF] Leach | 05 -Ligand-Based Approaches: Core Molecular Modeling[END_REF] 𝐸 𝑡𝑜𝑡𝑎𝑙 = 𝐸 𝑏𝑜𝑛𝑑𝑠𝑡𝑟𝑒𝑡𝑐ℎ𝑖𝑛𝑔 + 𝐸 𝑏𝑜𝑛𝑑𝑏𝑒𝑛𝑑𝑖𝑛𝑔 + 𝐸 𝑡𝑜𝑟𝑠𝑖𝑜𝑛𝑎𝑙 + 𝐸 𝑒𝑙𝑒𝑐𝑡𝑟𝑜𝑠𝑡𝑎𝑡𝑖𝑐 + 𝐸 𝑣𝑎𝑛𝑑𝑒𝑟𝑊𝑎𝑎𝑙𝑠 Eq. 58

Where, the first three terms (𝐸 𝑏𝑜𝑛𝑑𝑠𝑡𝑟𝑒𝑡𝑐ℎ𝑖𝑛𝑔 , 𝐸 𝑏𝑜𝑛𝑑𝑏𝑒𝑛𝑑𝑖𝑛𝑔 𝑎𝑛𝑑 𝐸 𝑡𝑜𝑟𝑠𝑖𝑜𝑛𝑎𝑙 ) represent the intramolecular terms and the two last terms (𝐸 𝑒𝑙𝑒𝑐𝑡𝑟𝑜𝑠𝑡𝑎𝑡𝑖𝑐 𝑎𝑛𝑑 𝐸 𝑣𝑎𝑛𝑑𝑒𝑟𝑊𝑎𝑎𝑙𝑠 ) are intermolecular terms. The forms of some of the commonly used potentials are presented below.

Molecular mechanics energy terms

Bond stretching term is often described by a Hook's law harmonic oscillator equation.

It is sometimes described by a Morse potential. [START_REF] Leach | 05 -Ligand-Based Approaches: Core Molecular Modeling[END_REF] 𝐸 𝑏𝑜𝑛𝑑s = ∑ 𝑘 𝑖 2 (𝑙 𝑖 -𝑙 𝑖,0 ) 2 𝑏𝑜𝑛𝑑𝑠 Eq. 59

Where 𝑙 𝑖 and 𝑙 𝑖,0 represent the length of the stretched and reference bonds, respectively. 𝑘 𝑖 are the bond constants.

Angle bending term is applied to all three atoms sequentially bonded. As with the bond term, a quadratic Hook's law equation is often used to describe the bond bending. [START_REF] Leach | 05 -Ligand-Based Approaches: Core Molecular Modeling[END_REF] 𝐸 𝑎𝑛𝑔𝑙𝑒𝑠 = ∑ 𝑘 𝑖 2 (𝜃 𝑖 -𝜃 𝑖,0 ) 2 𝑎𝑛𝑔𝑙𝑒𝑠 Eq. 60

Where 𝜃 𝑖 and 𝜃 𝑖,0 represent the size of the distorted and reference angles 𝑖, respectively. 𝑘 𝑖 are the force constants.

The Bond rotation (torsion) term considers four atoms sequentially bonded. Torsional potential term is most often described by a cosine series expression. 62

𝐸 𝑡𝑜𝑟𝑠𝑖𝑜𝑛𝑠 = ∑ 𝑉 𝑛 2 (1 + cos (𝑛𝜔 𝑖 -𝛾)) 𝑡𝑜𝑟𝑠𝑖𝑜𝑛 Eq. 61
Where 𝜔 𝑖 is the torsion angle 𝑖, 𝑉 𝑛 is a constant. 𝑛 is the multiplicity which corresponds to the number of the energy minima as the bond is rotated through 360°. 𝛾 is the phase factor which determines where the torsion potential passes through its minimum value.

P a g e | 48

The non-bonding intermolecular interactions consist of the van der Waals (vdW) forces and electrostatic interaction. The Coulomb's law term is widely used for the attraction or repulsion energies between point charges (e.g. atomic charge). 62

𝐸 𝑒𝑙𝑒𝑐𝑡𝑟𝑜𝑠𝑡𝑎𝑡𝑖𝑐 = ∑ 𝑁 𝑖=1 ∑ 𝑞 𝑖 𝑞 𝑗 4𝜋𝜀 0 𝑟 𝑖𝑗 𝑁 𝑖>𝑗 Eq. 62
Where 𝑞 𝑖 and 𝑞 𝑗 are the charges of the interacting atoms 𝑖 and 𝑗, respectively, and the distance between them is donated by 𝑟 𝑖j . 𝜀 0 is the vacuum permittivity.

The van der Walls (vdW) forces are usually described by so-called Lennard-Jones (6-12) potential. 62

𝐸 𝑣𝑑𝑊 = ∑ 𝑁 𝑖=1 ∑ (4𝜀 𝑖𝑗 [( 𝜎 𝑖𝑗 𝑟 𝑖𝑗 ) 12 -( 𝜎 𝑖𝑗 𝑟 𝑖𝑗 ) 6 ]) 𝑁 𝑖>𝑗 Eq. 63
where 𝑟 𝑖j is the distance between the centers of the non-bonded atoms 𝑖 and 𝑗. Terms 𝜀 𝑖𝑗 𝑎𝑛𝑑 𝜎 𝑖𝑗 are the Lennard-Jones well depth and collision diameter (the separation for which the energy is zero), respectively.

The cross and out-of-plane bending are other contributions terms, which can be included in a force field to ensure that the predicted geometry and energy are consistent with the experimental observations. Cross terms consider the correlation between stretch and bend components.

Electronic Polarisation is the phenomenon referring to the distortion of electron density of a molecule under the influence of an external field. The primary effect of the external field is to induce a dipole in the molecule. The polarisation is usually accounted for by calculating the interaction between induced dipoles. The induced dipole (𝛍) is calculated by multiplying the atomic polarizability (α) by the electric field (𝑬) present at that nucleus: 62 𝛍 = α𝑬 Eq. 64

Molecules in solution

Implicit and explicit solvent models

So far, we have omitted the solvent that plays a very important role in chemistry. The interactions between the solute-solvent molecules impact the general chemistry of the P a g e | 49

molecule being studied. These interactions can modify energy and stability. In addition, properties relating to energy (vibrational frequencies, spectrum, etc.) can be also changed.

Therefore, models that describe the interactions between the solute and the solvent are needed. Solvent models can be broadly categorised into two solvation models, implicit and explicit solvent models, which are distinguished by whether containing or not discrete solvent molecules. In explicit solvent models, solvent molecules are explicitly described in the calculations and the molecular details of each solvent molecule are treated. In these models, specific interactions between the solute and the solvent and between solvent molecules can be described and the individual interactions can be traced, this is especially interesting when interactions such as hydrogen bonds have to be described. The major drawback of the explicit models is that they are time-consuming and computationally expensive. In addition, the solvent molecules can adopt different conformations and a dynamic approach must then be used. For the calculation of the electronic properties, a way of overcome the huge problem of computational time is to adopt a QM/MM (Quantum Mechanics/Molecular Mechanics) approach in which the solvent molecules are described by molecular mechanic which greatly reduces the CPU time.

Implicit models represent the solvent molecules using a dielectric continuum with the Continuum implicit solvent approaches are highly popular due to their simplicity and low cost. On the other hand, the drawback of these approaches is that they do not account for specific solute-solvent interactions such as hydrogen bonding. These approaches have a major limitation due to the response time of the solvent (dynamic and inertial part) to changes of the solute charge distribution (solute electronic density). That is, the response time of the dynamic The non-equilibrium and equilibrium solvation states can be illustrated by the description of four steps corresponding to a fluorescence experience, starting from an initial solute-solvent GS equilibrium state, excitation and relaxation before returning back to the final solute-solvent GS equilibrium state, considering both solvent polarisation and solute geometry relaxations, as follows:

1)

In the GS optimisation step, the solute and solvent are in an equilibrium state.

The initial solute geometry relaxes towards the GS equilibrium geometry. The solvent (even the inertial part) is instantly relaxed and equilibrated towards the solute GS equilibrium geometry.

2)

In the vertical electronic transition step (absorption), the solute and solvent are in a non-equilibrium state. The solvent is partially (only the dynamic component, the solvent electrons) equilibrated with the new excited solute charge distribution.

3)

In the ES optimisation step, the solute and solvent are in an equilibrium state.

The ES solute geometry relaxes towards its ES equilibrium geometry. The solvent (even the inertial part) is instantly relaxed and equilibrated towards the solute ES equilibrium geometry.

P a g e | 52

4)

In the vertical electronic transition step (emission), the solute and solvent reach again a non-equilibrium state as in the absorption step. The solvent is partially equilibrated with the new ground solute charge distribution.

For excited state (ES) calculations using PCM models within a TD-DFT formalism, several modelling methods can be used to describe the solvent-solute response. The most commonly used methods are Linear-Response (LR) [START_REF] Cossi | Time-dependent density functional theory for molecules in liquid solutions[END_REF][START_REF] Cammi | Fast evaluation of geometries and properties of P a g e | 205 excited molecules in solution: a Tamm-Dancoff model with application to 4dimethylaminobenzonitrile[END_REF] and State-Specific (SS). [START_REF] Cossi | Solvent effect on vertical electronic transitions by the polarizable continuum model[END_REF] In order to perform properly the vertical electronic transition steps (non-equilibrium state), one needs to use the appropriate method, especially when a large change in a solute charge distribution is found. This is related to the difference between LR and SS formalisms. In the LR formalism, the response of the dynamic (electronic) part of the solvent is computed using the transition density while in the SS formalism, is computed as the difference between the electron densities of the final and the initial states. Moreover, along the vertical transition (non-equilibrium state), the SS formalism considers a mutual relaxation between the solute charge distribution and the dynamic part of the solvent response, which is ignored in the LR formalism. On the other hand, same inertial solvent response is included in both formalisms.

The LR method is generally more adequate to electronic transitions characterized by a small change of the solute charge distribution, since the electrostatic interactions between the solvent and the final solute charge distribution are not very different from those of the initial one. On the other hand, for electronic transitions involving large changes of the solute charge distribution, the LR method is insufficient because it does not consider the density-dependent relaxation of the solvent polarisation. In these cases, the excitation energies computed with the LR and the SS methods can be very different. Both methods provide the same final expressions for vertical transition energies when applied to isolated solutes.

In general, PCM introduces an effective solute Hamiltonian, in which an explicit solvent operator (known also as "reaction field" operator) is added to the solute Hamiltonian. Since the reaction field operator is a function of the solute charge distribution, nonlinear effects can then be induced in the effective Hamiltonian, which complicate the description of the transition energies. A new approach has been developed to overcome this intrinsic limit of the nonlinear effective Hamiltonian by introducing state specific solvent responses. It known as state-specific corrected linear response (cLR) approach. [START_REF] Caricato | Formation and relaxation of excited states in solution: A new time dependent polarizable continuum model based on time dependent density functional theory[END_REF] This cLR approach refines and corrects the linear response (LR) excitation energies.

P a g e | 53

Periodic boundary conditions

In order to set up a molecular dynamics simulation in explicit solvent it is needed to choose some appropriate boundary conditions to simulate processes in a small part (also called a unit cell or image of the simulation box) of a large system. [START_REF] Butler | An Introduction To Computational Methods[END_REF] When simulating a large system, a significant fraction would be close to the wall of the box. In order to avoid edge effects, boundary conditions are recommended. The most common used choice is periodic boundary conditions (PBC). In the PBC, each atom in the image of the simulation box (represented in grey in Figure 16) is replicated to form an infinite lattice filling the whole space. Actually, all atoms at the edge of the next box are a copy of the ones at the opposite edge of the box. Thus, all periodic images of the molecules will move exactly the same way the initial molecules in the central simulation box. In order to keep a fixed number of molecules in the central box, each leaving images will be systematically replaced by its image entering from the opposite face. On the other hand, since the box is replicated infinitely in all directions, van de Waals and electrostatic interactions exist not only between each nonbonded pair of atoms of the central box, but also with their mirror image in the neighbouring boxes. To overcome the problem of having an infinity of terms for these interactions, nonbonding interactions are spatially truncated using a cut-off (Rcutoff) distance of the interaction potential. In general, it required that the size of the simulation box should be larger than 2Rcutoff. In this way, each atom will have only interactions with those included inside the spherical cut-off (Figure 17), [START_REF] Haberlandt | Statistical Theory and Molecular Dynamics of Diffusion in Zeolites[END_REF] cannot interact with its mirror image, and will not be considered twice in the same interactions. This condition is called the minimum image criterion.

P a g e | 54

Water models

The water molecules can interact with the system of interest and form interactions such

Molecular mechanics force fields

The mainly difference between force fields resides in the number of terms in the energy expression, the complexity of those terms, and the way in which the constants were determined.

Commonly used molecular mechanics force fields are: AMBER [START_REF] Corongiu | [END_REF][74][75] (Assisted Model Building with Energy Refinement), CHARMM 76 (Chemistry at Harvard Macromolecular Mechanics), GROMOS 77,78 (Gronigen Molecular Simulation). The AMBER and CHARMM force fields were parameterized specifically for protein systems and nucleic acids. AMBER force field uses only five bonding and non-bonding terms along with a sophisticated electrostatic treatment, but the cross terms are neglected. The GROMOS force field is widely used for modelling biomolecules and for predicting the dynamical motion of molecules and bulk liquids. GROMOS and CHARMM force field uses five valence terms, one of which is an electrostatic term. The drawback of these three force fields is that they do not include the explicit polarisation effect, although it is important when simulating some systems. AMOEBA 79,80 (Atomic Multipole Optimised Energetics for Biomolecular Applications) force field is a polarisable force field, which uses atomic induced dipole to model polarisation and atomic multipoles up to quadrupoles to describe the permanent electrostatics. AMOEBA force field uses five terms to describe the short-range valence interactions (bond stretching, angle bending, torsional rotation, bond-angle cross term and out-of-plane bending), and three terms to describe the non-bonded vdW and electrostatic contributions. In this Ph.D. Thesis, we have used AMBER force field for modelling the protein system.

Energy minimisation

Energy minimisation is a local optimisation technique, used to bring the initial structure to the closest minimum on the potential energy surface (PES) by minimising its energy. In a potential energy landscape there generally exists more than one minimum, local and global minima, as represented in Figure 18. The global minimum being the minimum of the minima. [START_REF] Young | Computational Chemistry: A Practical Guide for Applying Techniques to Real-World Problems[END_REF][START_REF] Butler | An Introduction To Computational Methods[END_REF] Figure 18 1D schematic representation of the energy of all possible conformations of a simple molecule.

P a g e | 56

In order to find this global minimum of the energy hence, the most stable conformation, some algorithms need to be used to sufficiently sample the potential energy landscape. Most of these algorithms use the first and the second derivatives of the potential energy with respect to the geometric parameters. The most commonly used optimisation algorithms are the steepest descent, [START_REF] Arfken | The method of steepest descents[END_REF][START_REF] Petrova | The origin of the method of steepest descent[END_REF] the conjugate gradient [START_REF] Fletcher | Function minimization by conjugate gradients[END_REF] and the Newton-Raphson 84,85 algorithms.

The steepest descent algorithm

The steepest descent (also known as gradient descent) algorithm for energy minimisation considers only the first-order derivative information of the energy surface (the gradient 𝜵𝑓) to determine the search direction and the step size (or how far to descent). This means it only considers the first derivative when performing the updates on the geometric parameters. At each iteration step, the steepest descent method calculates the direction of the gradient. The steepest descent method can be written as

𝒙 𝑛+1 = 𝒙 𝑛 -𝛾𝜵𝑓(𝒙 𝑛 ) , 𝑛 ≥ 0 Eq. 66
where the step size (𝛾) is allowed to change at each iteration step to reach the local minimum.

This method starts at an arbitrary point (𝒙 0 ) and moves from 𝒙 𝑛 to 𝒙 𝑛+1 by minimising along the vector 𝒗 𝑛 extending from 𝒙 𝑛 in the direction of steepest descent, -𝜵𝑓(𝒙 𝑛 ), ensuring the steepest energy decrease. The vector (𝒗 𝑛 ) represents the search direction for minimising the objective function along the steepest descent direction, and it can be expressed as [START_REF] Chang | e-Design: Computer-Aided Engineering Design[END_REF] 𝒗 𝑛 = 𝒗(𝒙 𝑛 ) = -𝜵𝑓(𝒙 𝑛 ) ||𝜵𝑓(𝒙 𝑛 )|| Eq. 67

The steepest descent method takes a series of optimisation steps (𝒙 1 , 𝒙 2 , … , 𝒙 𝑛 ) until it is close enough to the minimum. [START_REF] Refsns | A brief introduction to the conjugate gradient method[END_REF] We consider that the minimum is reached when the difference in the objective function values of two consecutive iterations is less than a given convergence tolerance (threshold).

P a g e | 57

The steepest descent algorithm does not include the second-order derivative information. For this reason, the minimisation efficiency is not so good and the convergence rate is relatively slow close to the minimum. In addition, the convergence can fail for surfaces with long narrow valleys. However, the steepest descent algorithm is very simple and fast so, it is widely used to refine a structure that is far from the stationary point.

Conjugate gradient method

The conjugate gradient method is an improvement of steepest descent method. This method starts at an initial design (or starting point 𝒙 0 ), set the convergence tolerance, and calculate the function value (𝑓(𝒙 0 )) and gradient (𝜵𝑓(𝒙 0 )). [START_REF] Chang | e-Design: Computer-Aided Engineering Design[END_REF] Once the gradient is obtained, the search direction can be calculated for first iteration step. In general, the conjugate gradient method needs only few steps to find the minimum. Moreover, unlike the steepest descent method, the conjugate gradient method improves the convergence rate close to the minimum. [START_REF] Refsns | A brief introduction to the conjugate gradient method[END_REF] Figure 20 Illustration of minimisation energy using the conjugate gradient algorithm. Here, the minimum is reached from the starting point (𝒙 0 ) after two (𝒙 2 ) optimisation steps, with the steps of the conjugate gradient method in green.

Newton-Raphson method

The steepest descent method considers only the first-order derivative information to determine the search direction. In order to represent the objective function more accurately and to find a better search direction, hence to improve the optimisation efficiency, methods that include the second-order derivative information should be used. The inclusion of the second-order derivative information provides a better convergence rate as well. An example is the Newton-Raphson method, which uses the Hessian of the function in calculation the search direction. This method has a quadratic convergence rate, that is, it converges rapidly when the design point is within a certain radius of the minimum point.

P a g e | 58

Figure 21 Illustration for geometry optimisation using the Newton-Raphson.

In the Newton-Raphson method, at each optimisation step, the objective function is approximated by a quadratic function and moves towards its minimum upon reaching the minimum point (Figure 21). This method starts to find the second-order approximation of the objective function in the starting point (𝒙 0 ), then, find the minimal value location of it (𝒙 1 ). This can be determined by taking first differential, and make it equal to zero. By iterating this process, the minimum point can be reached.

It is noteworthy that saddle points can be determined using this method, which is infeasible using the steepest descent method. In addition, the number of the optimisation steps is reduced thanks to the evaluation of the second-order derivative information, unlike the steepest descent method. The main drawback of this method lies to the step size, which goes towards infinity when the curvature of the potential energy surface is close to zero. Methods using both the gradient and second-order derivative (Hessian matrix) often require fewer minimisation steps but more CPU time due to the time needed to compute the Hessian matrix.

Molecular dynamics (MD)

Molecular dynamics (MD) is a computational method, which simulates the dynamic behaviour of a molecular systems. It can be used to study conformational space, ligand binding to enzymes, small molecules adsorbing to a surface, and is often the method of choice for large systems. Molecular dynamics trajectories (series of atomic positions in time) of atoms are determined by numerically solving Newton's laws of motion (classical equations of motion) for a system of interacting particles, where forces between the particles and potential energy are represented by molecular mechanics force fields.

Having the force (𝐹 𝑖 (𝑡)) on each atom (𝑖) at time 𝑡, the acceleration (𝑎 𝑖 ) can be computed from the force acting on atom 𝑖 according to Newton's second law, and forces can be calculated from the molecular mechanics force field, according to the following equations: [START_REF] Butler | An Introduction To Computational Methods[END_REF] P a g e | 59

𝐹 𝑖 (𝑡) = 𝑚 𝑖 𝑎 𝑖 (𝑡) Eq. 68

where 𝑚 𝑖 is the mass of the atom. This way, the integration of the motion equations over time provides the new particle positions, velocities and accelerations at each time moment (time step).

𝐹 𝑖 (𝑡) = -𝑚 𝑖 𝑑𝑉(𝑟(𝑡)) 𝑑𝑟 𝑖 (𝑡) = 𝑚 𝑖 𝑑 2 𝑟 𝑖 (𝑡) 𝑑𝑡 2
Eq. 69

where 𝑟 𝑖 (𝑡) represents the vector position of each atom 𝑖 at time 𝑡. 𝑉(𝑟(𝑡)) is the potential energy function of the system. In particular, time is partitioned into time steps (𝛿𝑡), which are employed to propagate the system forward in time. In order to determine a trajectory of each atom it is needed to know the initial position of the atoms, an initial distribution of velocities and the acceleration. Diverse integration algorithms have been developed to integrate Newton's equations of motion by discrete-time numerical approximations. The commonly used algorithms to propagate coordinates over time are the Verlet, [START_REF] Verlet | [END_REF] Verlet Leapfrog [START_REF] Rw | Computer simulation using particles[END_REF][START_REF] Van Gunsteren | A leap-frog algorithm for stochastic dynamics[END_REF] and Velocity-Verlet [START_REF] Wc | A computer simulation method for P a g e | 206 the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters[END_REF][START_REF] Mp. Allen | Computer simulation in chemical physics[END_REF] algorithms. These algorithms assume that the positions, velocities and accelerations can be approximated as a second order Taylor expansion.

Verlet algorithm

The Verlet algorithm takes the positions and accelerations of the atoms at the current time step and positions from the previous step to calculate the positions for the next time step:

𝑟 𝑖 (𝑡 + 𝛿𝑡) = 2𝑟 𝑖 (𝑡) -𝑟 𝑖 (𝑡 -𝛿𝑡) + 𝑎 𝑖 (𝑡)𝛿𝑡 2 Eq. 70

where 𝑟 𝑖 (𝑡) and 𝑎 𝑖 (𝑡) are the position and acceleration of atom 𝑖 at time 𝑡, respectively. 𝛿𝑡 is the time step. The 𝑟 𝑖 (𝑡 -𝛿𝑡) and 𝑟 𝑖 (𝑡 + 𝛿𝑡) are the previous and next positions for given time step, respectively. The Verlet algorithm assumes that the forces depend uniquely on the positions. In this algorithm, the velocity is not treated explicitly, however, it is needed to estimate the kinetic energy.

LeapFrog algorithm

The LeapFrog algorithm is obtained by a simple modification of the Verlet integrator. In this Leapfrog algorithm, the positions (𝑟) at a time (𝑡 + 𝛿𝑡) are computed through the velocities at time 𝑡 + where 𝑟 𝑖 (𝑡) and 𝑎 𝑖 (𝑡) are the position and acceleration of atom 𝑖 at time 𝑡, respectively. 𝑣 𝑖 is the velocity of atom 𝑖.

The velocity is not treated accurately in both Verlet and LeapFrog integration schemes.

In order to handle the velocities in satisfactorily manner, velocity-Verlet algorithm has been developed.

Velocity-Verlet algorithm

This algorithm takes positions, velocities, and accelerations at the current time step. This gives a more accurate integration than the Verlet algorithm. Velocity-Verlet integrator computes position and velocity of an atom (𝑖) at the time step (𝑡 + 𝛿𝑡), starting from step (𝑡):

𝑟 𝑖 (𝑡 + 𝛿𝑡) =

Thermodynamic ensembles

Molecular dynamics simulations can be performed in various thermodynamic

ensembles. An ensemble is assumed as an imaginary collection of systems described by the common microscopic interactions and sharing a common set of macroscopic properties.

Thermodynamic ensemble is distinguished by which thermodynamic variables are kept constant over the course of the simulation. Among the thermodynamic ensembles used in MD simulations is the Microcanonical Ensemble (NVE). This ensemble keeps constant the number of particles (N), the volume of the simulation cell (V), and the total energy of the system (E). This is defined as a default in MD simulation. Other relevant ensembles are the Canonical Ensemble (NVT) which describes a system at constant particles number (N), the volume of the simulation cell (V), and temperature (T) and the Isothermal-Isobaric Ensemble (NPT) in which number of particles (N), temperature (T), and pressure (P) are kept constant, while the volume and the energy are allowed to fluctuate. A constant temperature and pressure are maintained by using thermostat and barostat. The most rigorous methods to control the temperature are Berendsen, [START_REF] Berendsen | Molecular dynamics with coupling to an external bath[END_REF] Andersen 93 and Nosé-Hoover [START_REF] Nosé | A unified formulation of the constant temperature molecular dynamics methods[END_REF][START_REF] Hoover | Canonical dynamics: Equilibrium phase-space distributions[END_REF] thermostats. In NPT ensemble, the temperature is controlled using Langevin dynamics in which the Berendsen thermostat is accompanied with a barostat, whereby a piston controls the fluctuations of the volume and hence maintains constant pressure. In this Ph.D. Thesis, the MD simulations (using periodic boundary conditions)

have been performed under NPT ensemble (T = 300 K and P = 1 atm). During these simulations, the pressure and temperature have been maintained by the Berendsen algorithm. [START_REF] Berendsen | Molecular dynamics with coupling to an external bath[END_REF]

Umbrella Sampling (US)

Free energy is the measure for the driving force of any process, such as a chemical reaction. Among the approaches designed to calculate the free energy is the Umbrella Sampling (US) method. [START_REF] Torrie | Nonphysical sampling distributions in Monte Carlo free energy estimation: umbrella sampling[END_REF] The US method is often used to compute the potential of mean force (PMF) of noncovalent ligand-receptor association/dissociation systems or to compare the stability of two different conformations and identify transition states. The advantage of this method is its ability to sample a large number of conformations, which some of them are inaccessible during a classical MD. The basic idea of the US method is to introduce the concept of biased simulation (usually called window) in order to generate an enhanced sampling over a focused region of configuration space, which is achieved by applying an external biasing potential (restraint) for each window (for each simulation). This restraint forces the reaction coordinate to remain close to the center of the window. The samplings in each window should have some overlap with P a g e | 62

adjacent windows. The information from the different biased windows is converted into a series of histograms, which contain the biased distribution of the reaction coordinates. These histograms are then unbiased and combined, usually using the weighted histogram analysis method [START_REF] Kumar | The Weighted Histogram Analysis Method for Free-Energy Calculations on Biomolecules.1, The Method[END_REF][START_REF] Roux | The Calculation of the Potential of Mean Force Using Computer-Simulations[END_REF] (WHAM), to generate the Potential of Mean Force (PMF).

Molecular Docking

Docking 99,100 is a molecular modelling approach widely used for the study of molecular recognition, which aims to identify correct conformation/orientation of ligands inside a protein binding pocket and to predict the binding affinity between the protein and the ligands. Docking process permits to estimate the binding energy, free energy and the stability of complexes.

According to the type of ligand, the docking can be classified as protein-ligand docking, proteinprotein docking, or protein-nuclei acid docking. For performing molecular docking, a 3D of the molecular target is one of the necessary conditions, which can be an experimentally solved structure (such as NMR and X-ray crystallography structures) or homology-modelled structure (obtained by computational techniques). Docking can be achieved through two steps: i) an engine for conformations sampling of the ligand in the protein active site and ii) scoring function, which associates a score to each predicted conformation and then correctly ranks these conformations.

P a g e | 63

Quantum mechanics/molecular mechanics (QM/MM) method

The hybrid quantum mechanics/molecular mechanics (QM/MM) method has become a popular method for simulating and investigating reactions inside a large system such as proteins (enzymes). The idea behind the hybrid QM/MM method was introduced by Warshel and Levitt 101

to simulate enzymatic reactions. The QM/MM method was designed by combining both quantum mechanics (QM) and molecular mechanics (MM) methods into one calculation. The QM method can model chemical reactions and compute many properties (e.g., electronic excitations, charge transfer), but it is computationally infeasible for large systems. The MM method is able to treat and describe large systems (several thousand particles) using molecular mechanics force fields. However, these force fields do not treat the electronic part of the system.

Therefore, the entire system is partitioned into two subsystems, and treats with accurate and computationally expensive methods only the necessary region (the "QM subsystem"), while the rest is described at lower level of theory (the "MM subsystem") (Figure 24). The QM subsystem corresponds to the active site of the protein in where the phenomenon of chemical interest is taking place. This QM subsystem is treated quantum-mechanically. MM subsystem corresponds to the reminder of the system, the surrounding environment of the QM subsystem. This MM subsystem is described by MM force fields.

P a g e | 64

Figure 24 illustration of the QM/MM system. The MM subsystem represents the environment (in red: the water molecules, in ribbon: the protein). The QM subsystem represents the substrate.

QM/MM method is widely used as a valuable tool for the study of explicit solvent and biomolecular processes. In this Ph.D. Thesis, the QM/MM method has been applied as follow:

the substrate (oxyluciferin) present inside the active site of the protein (luciferase) is treated at the QM level while the remainder of the system (protein, solvents…) is described by parameters listed in the MM force fields. In this hybrid QM/MM method, the total potential energy of the system can be expressed as:

𝐸 𝑡𝑜𝑡𝑎𝑙 𝑄𝑀/𝑀𝑀 = 𝐸 𝑄𝑀 + 𝐸 𝑀𝑀 + 𝐸 𝑄𝑀-𝑀𝑀 Eq. 75
where 𝐸 𝑄𝑀 represents the interactions between the atoms of the QM subsystem, 𝐸 𝑀𝑀 represents the interactions between the atoms of the MM subsystem, and 𝐸 𝑄𝑀-𝑀𝑀 describes the interactions between the atoms of the QM and MM subsystems. This last term is more difficult to describe, and in order to overcome this difficulty several QM/MM schemes have been developed, basically differing on the form of the 𝐸 𝑄𝑀-𝑀𝑀 . These QM/MM schemes can be divided into two distinct categories: subtractive and additive coupling schemes.

Subtractive and additive QM/MM coupling schemes

Subtractive QM/MM coupling schemes

In the subtractive QM/MM coupling scheme, the QM/MM energy of the system is determined in three steps. The first step evaluates the total energy of the entire system, consisting of both QM and MM subsystems, at the MM level. The second step determines the QM energy of the QM subsystem at the QM level of theory. The third step computes and subtracts the MM energy of the QM subsystem in order to ensure that the contribution of the QM subsystem energy does not included twice to the total energy as it has already been computed in the first step at the MM level of theory. The advantage of this QM/MM scheme is that no explicit coupling terms are required between the two subsystems, which make the implementation fairly straightforward. On the other hand, MM force field is required for the QM subsystem, which is not generally available and can be difficult to obtain. In the original ONIOM method, the coupling terms between the two subsystems are described entirely at the MM level, this makes the subtractive QM/MM scheme problematic for the electrostatic interactions, which are described by fixed atomic charge at the MM level. 102 Several improvements of the original ONIOM method were recently made, allowing to include the MM charges in the QM Hamiltonian. 102,103

Additive QM/MM coupling schemes

In the additive QM/MM coupling scheme, the QM subsystem is embedded within the larger MM subsystem. The QM/MM energy of the system is determined by calculating not only the sum of the energy terms of MM and QM subsystems but also the coupling energy terms between the two subsystems (interactions energy between the substrate and the surrounding).

The energy expression for the additive QM/MM scheme can be written as [START_REF] Monticelli | Biomolecular Simulations: Methods and Protocols[END_REF] 𝐸 𝑡𝑜𝑡𝑎𝑙 𝑄𝑀/𝑀𝑀 = 𝐸 𝑀𝑀 𝑀𝑀 𝑙𝑒𝑣𝑒𝑙 + 𝐸 𝑄𝑀 𝑄𝑀 𝑙𝑒𝑣𝑒𝑙 + 𝐸 𝑄𝑀+𝑀𝑀 𝑐𝑜𝑢𝑝𝑙𝑖𝑛𝑔 𝑄𝑀-𝑀𝑀 Eq. 77

In this additive QM/MM scheme, the interactions between the two subsystems are described explicitly by a coupling term (𝐸 𝑄𝑀+𝑀𝑀 𝑐𝑜𝑢𝑝𝑙𝑖𝑛𝑔 𝑄𝑀-𝑀𝑀 ). Contrary to the subtractive QM/MM scheme, only the interactions within the MM subsystem are described at the MM level, giving the 𝐸 𝑀𝑀 𝑀𝑀 𝑙𝑒𝑣𝑒𝑙 while the QM subsystem is treated at the QM level as before (𝐸 𝑄𝑀 𝑄𝑀 𝑙𝑒𝑣𝑒𝑙 ).

Mechanical embedding

In a mechanical embedding model, the QM-MM interactions between the subsystems are treated as the MM-MM interactions, that is, at the MM level. 102,104,105 The QM subsystem is thus kept in place by MM interactions. These interactions include both bonded (bond stretching, bond bending, and torsional) interactions and non-bonded (electrostatic and van der Waals) interactions. In the simplest implementation of mechanical embedding model, the electronic wavefunction is evaluated for an isolated QM subsystem. Therefore, the MM subsystem cannot polarise the electron density in the QM subsystem. In order to compute the electrostatic interactions between the QM and MM subsystems, one can use a fixed set of charges for the QM subsystem obtained by the force field. The mechanical embedding model is the lowest level of interaction as it does not include the polarisation and only bonded and steric energies are considered.

Electrostatic embedding

Electrostatic embedding is an improvement of the mechanical embedding model in where 𝑞 𝐵 represents the MM point charge and 𝑍 𝐴 is the nuclear charge of the QM atoms. 𝑁 is the number of electrons. 𝑖, 𝐵 and 𝐴 run over all the electrons, point charges and QM nuclei, respectively.

In such an embedding model the electronic structure of the QM subsystem can adapt to changes in the charge distribution of the environment and is systematically polarised by it. The QM-MM electrostatic interactions are treated at the QM level, which gives a more advanced and accurate description than a mechanical embedding model. Hence, this can lead to a high computational cost, especially for the computation of the electrostatic force due to the QM density acting on the MM point charges. Special care is required at the boundary between the two subsystems, where the MM point charges are very close to the QM electron density, which can cause an over-polarisation. This problem is particularly pronounced when the QM-MM boundary crosses a covalent bond. In this Ph.D. Thesis, we have applied the Electrostatic embedding.

Polarisation embedding

In an electrostatic embedding, the interaction of the polarisable QM density with MM point charges is considered. Now, the next step is to include the polarisation of the MM atoms by the QM charges distribution. In the polarisation embedding model both QM and MM subsystems can mutually polarise each other. 105 Thus, not only the atoms of the MM subsystem can polarise the QM subsystem, but also the MM subsystem can be polarised by the atoms of the QM subsystem. In the polarisation embedding model the total QM/MM energy is obtained by calculating the polarisation of the MM atoms at every step of the self consistent-field iteration of the QM wavefunction. The use of this embedding model offers a more realistic coupling between the subsystems. On the other hand, the implementation of polarisation force fields for biomolecular system is not straightforward. subsystem. In order to overcome this problem several approaches have been proposed. The most commonly used are:

Covalent bonds across the QM/MM boundary

Link atom approach

This approach consists in introducing an additional monovalent link atom at an appropriate position along the bond vector between the atoms of the QM and MM subsystems (Figure 27). For instance, the most often used atom to cap the QM subsystem is the hydrogen atom. [START_REF] Monticelli | Biomolecular Simulations: Methods and Protocols[END_REF]102 The monovalent link atom is present only in the QM calculations, and is invisible for the atoms of the MM subsystem.

Localised orbitals approach

This approach consists in replacing the chemical bond at the QM-MM boundary by a double-occupied molecular orbital. Several methods have been elaborated to introduce these orbitals. One of the most used method is: the local self-consistent field (LSCF) method, which introduces orbitals at the QM atoms. In the LSCF method, the atomic orbitals on the QM atom of the broken bond are localised and hybridised, while the hybrid orbital pointing to the MM atom is double-occupied (Figure 28). The other orbitals are single-occupied. During the SCF iterations of the QM wavefunction, the double-occupied molecular orbitals are kept frozen (i.e., they do not participate in the SCF iterations), while the other hybrid orbitals are optimised along with all orbitals in the QM subsystem. Chapter 4

Vibrationally-resolved spectra

In the following part, the procedure that permit to generate vibrationally-resolved electronic spectrum will be reviewed. We start by introducing the transition dipole moment integral followed by representing the Franck-Condon (FC) and Herzberg-Teller (HT) approximations.

The electronic spectrum such as ultraviolet-visible (UV-VIS) spectrum is simulated by computing vertical transition (excitation or de-excitation) energies at optimised geometries for the ground or excited electronic states. The resulting spectrum is usually simulated as a Gaussian convolution of the excitation (or de-excitation) energies, using a half-width at half-maximum (HWHM). Line intensities in one-photon absorption or emission spectra mainly depends on the square of the transition dipole moment (〈𝜓 ′ | 𝜇 | 𝜓 ′′ 〉) (1) and on the radiation frequency 𝜔:

𝜎 𝑎𝑏𝑠 (𝜔) = 4𝜋 2 𝜔 3𝑐 ∑ 𝜌 ′ |〈𝜓 ′ | 𝜇 | 𝜓 ′′ 〉| 2 ′,′′
𝛿(ℏ𝜔 -∆𝐸 -𝐸 ′′ + 𝐸 ′ ) Eq. 79 (1) The single and double quotes (′ 𝑎𝑛𝑑 ′′) indicate the initial and final states, respectively. For the absorption, the initial (single quotes) and final states (double quotes) represent the ground and the excited states, respectively. For the emission, the initial (single quotes) and final states (double quotes) represent the excited and the ground states, respectively.

𝜎 𝑒𝑚𝑖 (𝜔) = 4𝜔 3 3𝑐 3 ∑ 𝜌 ′ |〈𝜓 ′ | 𝜇 | 𝜓 ′′ 〉| 2 ′,′′
𝛿(∆𝐸 + 𝐸 ′ -𝐸 ′′ -ℏ𝜔) Eq. 80

Where 𝜎(𝜔) is the rate of photon absorption (or emission) per a single molecule and per unit of radiant energy. 𝑐 is the velocity of light and 𝜌 ′ is the Boltzmann population of a given initial state. 𝜓 ′ and 𝜓 ′′ are the molecular wavefunctions of the initial and final states, respectively. 𝜇 is the electric dipole moment. ∆𝐸 represents the energy difference between the minima of the potential energy surfaces (PESs). 𝐸 ′ or 𝐸 ′′ is the vibrational energy in the corresponding electronic state.

Transition dipole moment integral

According to the Born-Oppenheimer approximation, the molecular wavefunction of each state can be written as a product of an electronic (𝜓 𝑒 ) and a nuclear (𝜓 𝑛 ) wavefunctions.

Therefore, the transition dipole moment can be expressed in two terms, electronic (𝜇 𝑒 ) and nuclear (𝜇 𝑛 ) terms, as:

〈𝜓 ′ | 𝜇 | 𝜓 ′′ 〉 = 〈 𝜓 𝑒 ′ 𝜓 𝑛 ′ | 𝜇 𝑒 | 𝜓 𝑒 ′′ 𝜓 𝑛 ′′ 〉 + 〈𝜓 𝑒 ′ 𝜓 𝑛 ′ | 𝜇 𝑛 |𝜓 𝑒 ′′ 𝜓 𝑛 ′′ 〉 Eq. 81
It should be noted that the nuclear term (𝜇 𝑛 ) is independent from the electrons coordinates and hence, the second term of the previous equation can be expressed as:

〈𝜓 𝑒 ′ 𝜓 𝑛 ′ | 𝜇 𝑛 | 𝜓 𝑒 ′′ 𝜓 𝑛 ′′ 〉 = 〈𝜓 𝑛 ′ | 𝜇 𝑛 | 𝜓 𝑛 ′′ 〉 〈 𝜓 𝑒 ′ | 𝜓 𝑒 ′′ 〉 Eq. 82
Due to the orthogonality of the electronic wavefunctions to each other, terms of the different electronic states will be zero (〈 𝜓 𝑒 ′ | 𝜓 𝑒 ′′ 〉 = 0). Consequently, the second term in the equation (Eq. 81) will vanish. The transition dipole moment (Eq. 81) can then be expressed as:

〈𝜓 ′ | 𝜇 | 𝜓 ′′ 〉 = 〈 𝜓 𝑒 ′ 𝜓 𝑛 ′ | 𝜇 𝑒 | 𝜓 𝑒 ′′ 𝜓 𝑛 ′′ 〉 Eq. 83
or in more compact formula

〈𝜓 ′ | 𝜇 | 𝜓 ′′ 〉 = 〈 𝜓 𝑛 ′ | 𝜇 𝑖𝑓 | 𝜓 𝑛 ′′ 〉 with 𝜇 𝑖𝑓 = 〈 𝜓 𝑒 ′ | 𝜇 𝑒 | 𝜓 𝑒 ′′ 〉 Eq. 84
Now, the transition dipole moment depends on the nuclear wavefunctions (𝜓 𝑛 ′ , 𝜓 𝑛 ′′ ) and on the electronic transition moment (𝜇 𝑖𝑓 ). Within Eckart conditions, 112 the nuclear wavefunction can be written as: 𝜓 𝑛 = 𝜓 𝑟 𝜓 𝑣 where 𝜓 𝑟 and 𝜓 𝑣 are the rotational and the vibrational wavefunctions, respectively. If the rotation energies as sufficiently small, one can express the transition dipole moment as:

〈𝜓 ′ | 𝜇 | 𝜓 ′′ 〉 ≃ 〈 𝜓 𝑣 ′ | 𝜇 𝑖𝑓 | 𝜓 𝑣 ′′ 〉.

Franck-Condon (FC) approximation

Since the analytic expression of the electronic transition dipole moment (𝜇 𝑖𝑓 ) is remained unknown, an appropriate approximation is needed to compute the transition dipole moment. An approximation has been proposed by Franck and Condon, known as Franck-Condon (FC) principle 113,114 (Figure 30). The FC principle supposes that during an electronic transition, the electron jump occurs in such a short time that the positions and velocities of the nuclei are mostly unaltered by the molecular vibrations. That is, the motions of the nuclei are ignored during the electronic transition. However, the assumption that the electronic transition dipole moment is constant during the electronic transition makes the FC approximation very restrictive.

Herzberg-Teller (HT) approximation

To overcome the limitation of FC approximation, another approximation relying on the FC principle has been adopted, in which a small variation of the relative positions and velocities of the nuclei during the electronic transition is considered. This approximation has been improved by expanding the electronic transition dipole moment in a Taylor series of the normal coordinates about the equilibrium geometry of either the initial or the final electronic state. For instance, The Taylor expansion on the equilibrium geometry of the final electronic state can be written as: Here the Taylor development is expanded up only to the second derivatives. The zeroth-order term refers to a static electronic transition dipole moment. It is referred as the FC approximation. The first-order term considers a small variation of the dipole moment during the electronic transition. It is referred as the Herzberg-Teller (HT) 115,116 approximation (Figure 30).

𝜇 𝑖𝑓 (
𝑄 𝑒𝑞 ′′ represents the equilibrium geometry of the final state. N is the number of normal modes. If the fully-allowed electronic transitions (𝜇 𝑖𝑓 (𝑄 𝑒𝑞 ′′ ) ≫ 0) is reached, the overall spectrum will be described correctly. If the transition is not allowed or forbidden (𝜇 𝑖𝑓 (𝑄 𝑒𝑞 ′′ ) ≃ 0), the most intense vibronic transitions will be missed and hence, the overall spectrum will not be described properly. 𝑄 ′′ is the normal coordinates of the final state, which represents the all nuclear motions for a given vibration.

P a g e | 72 From the equation (Eq. 85), the transition dipole moment integral can then be expressed as:

〈 𝜓 ′ | 𝜇 | 𝜓 ′′ 〉 = 𝜇 𝑖𝑓 (𝑄 𝑒𝑞 ′′ )〈 𝜓 𝑣 ′ | 𝜓 𝑣 ′′ 〉 + ∑ ( 𝜕𝜇 𝑖𝑓 𝜕𝑄 𝑘 ′′ ) 0 〈 𝜓 𝑣 ′ |𝑄 𝑘 ′′ | 𝜓 𝑣 ′′ 〉 𝑁 𝑘=1 + 1 2 ∑ 𝑁 𝑘=1 ∑ ( 𝜕 2 𝜇 𝑖𝑓 𝜕𝑄 𝑘 ′′ 𝜕𝑄 𝑙 ′′ ) 0 〈 𝜓 𝑣 ′ |𝑄 𝑘 ′′ 𝑄 𝑙 ′′ | 𝜓 𝑣 ′′ 〉 𝑁 𝑙=1
Eq. 86

Where 〈 𝜓 𝑣 ′ | 𝜓 𝑣 ′′ 〉 is the overlap integral and it refers to the FC integral.

However, a problem will arise when trying to compute the overlap integral because each vibrational wavefunction can be expressed in a different set of normal coordinates. To overcome this, a solution has been proposed by Duschinsky,117 which expresses the normal coordinates of the initial state as a linear combination of the normal coordinates of the final state as:

𝑄 ′ = 𝐽𝑄 ′′ + 𝐾 Eq. 87

P a g e | 73
Here, 𝑄 ′ and 𝑄 ′′ are the normal coordinates of the initial and final states, respectively. The Duschinsky matrix 𝐽 represents the mixing of the normal modes during the transition. 𝐾 represents the displacement vector due to the changes in geometry structure between the initial and final states. Duschinsky matrix 𝐽 and displacement vector 𝐾 can be defined with the following relations:

𝐽 = (𝐿 ′ ) 𝑇 𝐿 ′′ , 𝐾 = (𝐿 ′ ) 𝑇 𝑀 1 2 𝑆 Eq. 88
where 𝐿 ′ and 𝐿 ′′ represent the transformation matrices from the Cartesian coordinates to normal coordinates, 𝑆 is the displacement vector of the nuclear Cartesian coordinates between the initial and final states, and 𝑀 represents the diagonal matrix of atomic masses.

In general, Duschinsky transformation is a good approximation, especially when the geometry structure does not undergo a major distortion during the transition.

The vibronic spectrum can be simulated within time-independent (TI) 118 or timedependent (TD) 119,120 frameworks. In the TI framework, the overall spectrum can be generated as an ensemble of distinct transitions between vibrational states of the initial and final electronic states. This is the framework of choice when high-resolution spectrum is available. [START_REF] Bloino | Aiming at an accurate prediction of vibrational and electronic spectra for mediumto-large molecules: an overview[END_REF] It suffers from convergence limitation related to the line shapes of absorption and emission spectra, especially when temperature effects are considered. To overcome this limitation, one needs to use the TD framework since it provides fully converged line shapes of absorption and emission spectra at finite temperature accounting for both Duschinsky and HT effects.

Since the initial and final state harmonic potential energy surfaces (PESs) can differ from each other and their minima can shift, the definition of the displacement vector (𝐾) (Eq. 87)

will become very challenging. There exist two models to describe the PES of the initial and final states, called adiabatic and vertical models. In adiabatic model, the PES of the initial and final states are constructed starting from a harmonic analysis at the corresponding equilibrium structures. This gives a best description of the equilibrium structures of the initial and final states and the frequencies. The vertical model emphasises the accuracy of the preeminent bands in the overall line shapes, thus the most intense transitions. However, it requires a better representation of the PES minimum about the vertical transition. In this model, the force constants (hessian matrix) are calculated in each electronic state.

Furthermore, the adiabatic and vertical models refer to as Adiabatic Hessian (AH) and Vertical Gradient (VG), respectively. AH model adopts the same electronic state harmonic PES both for absorption and emission and hence, is more adequate to analyse the difference of their

Generation of vibronic spectra

The FC, HT and FCHT approximations of the electronic transition dipole moment have been implemented in GAUSSIAN package. The vibrationally-resolved electronic spectrum can be generated by specifying a single keyword (FC or HT or FCHT) representing the terms of Taylor series (Eq. 85) which must be used in the calculation. An example of a simulated vibronic spectra (correspond to a one-photon absorption transition) using FC, HT and FCHT approximations with AH model is represented in Figure 31. The vibronic spectra can be simulated within TI or TD frameworks depending on whether or not temperature effects are included in the simulation. Furthermore, different models can be applied to describe the electronic transition such as AH, VG models. Each model recovers the frequencies and normal modes from the initial state. Depending on the features of the chosen model, the normal modes and frequencies from the final state can be ignored (e.g., VG model)

or considered (e.g., AH model).

In order to generate the vibrationally-resolved spectrum, one needs to prepare correctly the data files of the initial and the final states (equilibrium geometry, frequencies, normal modes…) and to specify keywords of the adequate approach. Whatever the chosen approach,

Gaussian package extracts the information from the data files. In addition, it maximises the overall superposition of the initial and final state structures in order to reduce the translation P a g e | 75

and rotation effects. The translation can be removed by superposing the center of masses of the two equilibrium geometries. The rotation is then minimised by maximising the superposition.

This procedure avoids numerical problems when the initial and final state geometries are slightly different. In a case where the difference between the two geometries is important (Figure 32), the superposition cannot be done correctly and hence, the vibronic spectrum simulation will not be generated when using some model such as AH model. It should be noted that the simulation of vibronic spectrum of superimposition geometries using the Adiabatic Hessian (AH) model can be done by forcing the calculations (by using specific commands in Gaussian programme). Therefore, the final AH vibronic spectrum cannot be much reliable as shown in the vibrationally-resolved absorption spectrum of protonated phenol-cycle where the electronic ground and first excited states minimum energy structures are superimposed (Figure 33a). On the other hand, for the deprotonated phenolcycle, both electronic ground and first excited states minimum energy structures are superposed and thus, AH method generates correctly the vibronic spectrum (Figure 33b).

P a g e | 76

If the superposition of the geometries is done successfully (natural oxyluciferin in Figure 32 and deprotonated phenol-cycle in Figure 33b), the relative transition energy with respect to the transition energy between the vibrational ground states of the initial and the final states (〈 0 ′ | 0 ′′ 〉), the square of the transition dipole moment integral and the intensity of the corresponding spectral line will be obtained. Furthermore, the broadening of the spectral line will be done by Gaussian functions with a certain HWHM value.

In case where the relative position of the ground state geometry of a given compound is not the same in the excited state, it is possible to apply VG model to generate vibronic spectrum. In order to get an idea on the involved effects when applying AH or VG models, we have performed vibronic calculations of two oxyluciferin analogues, phenol-OMe and OMe-enol using FCHT approximation with AH and VG models (Figure 34). The spectra shapes obtained with VG and AH models are very different. The spectra obtained with AH model represent a larger and more intense amount of vibronic transitions in the spectral region around 390-450 nm (Figure 34). In the vibronic spectra obtained with VG model these vibronic transitions (around 410-430 nm) are less intense. Furthermore, in both AH and VG models, the predicted 0-0' vibronic transition is the most intense one, except for the predicted one of phenol-OMe obtained with AH model. However, it is still close to the wavelength maximum of vibronic spectrum. Moreover, the 0-0' bands predicted with VG model are blue-shifted around 0.05 eV comparing with those obtained with AH model. This difference in energy is rather small, notwithstanding the VG model weaknesses. P a g e | 77

Conclusion

In this Chapter, we have presented a brief overview of the applied approaches to approximate the electronic transition dipole moment and generate the vibronic spectra. In addition, we have generated the vibronic spectra on the basis of the TD-DFT vibrational signatures by using the Franck-Condon (FC) and the Franck-Condon-Herzberg-Teller (FCHT) approximations within time independent (TI) framework and Adiabatic Hessian (AH) model as implemented in Gaussian 09. 109,111,121,122 We observed that the vibronic spectra obtained using FCHT approximation resemble the FC ones. We also simulated the vibronic spectrum using FCHT within TI framework and Vertical Gradient (VG) and compared it with the one simulated using Chapter 5

Computational details

In this Chapter, the computational procedures that have been used in this work will be presented in detail. This Chapter is sectioned as follows: We start by detailing the steps of the classical MD molecular dynamics performed for both compound in water and compound in protein, followed by describing the procedure that has been done to extract a statistical number of snapshots (conformations) for the QM/MM calculations to simulate the absorption and emission spectra. Moreover, the methods used for optimising the geometries and computing the absorption and emission transition (vertical transitions) energies are presented in the QM and QM/MM calculation sections. In addition, the computation details of the vibrationallyresolved spectra are reported in the vibronic spectra section.

Classical molecular dynamics (MD) simulations: compound in explicit water

A Classical molecular dynamics (MD) simulation of 10 ns has been performed with Amber14 program 123 for each system under study in order to obtain a statistical number of snapshots for performing QM/MM calculations to simulate the absorption (or emission) spectra considering explicit water molecules. In order to achieve this, the optimised geometry of the compound in the ground (or excited) state has been selected as the starting point for performing MD simulations. The bond distances, angles, dihedrals, and Mulliken atomic charges parameters for the compounds (at both ground and first excited states minima) have been taken from the QM calculations in gas phase and chosen as the starting point.

Each compound has been solvated using TIP3P water molecules model 124 within truncated octahedral box, ensuring a solvation shell of 10 Å around the compound, corresponding to around 800 solvent molecules.

Moreover, after the minimisation and heating (from 100 to 300 K in 20 ps) steps, a production step of 10 ns MD simulation (using periodic boundary conditions) has been P a g e | 80 performed with a 2 fs time step, under NPT (isothermal-isobaric ensemble) conditions (T = 300 K and P = 1 atm). During this simulation, the pressure and temperature have been maintained by the Berendsen algorithm. [START_REF] Berendsen | Molecular dynamics with coupling to an external bath[END_REF] Then, a statistical number of snapshots (100 conformations) have been extracted from the MD simulation for further QM/MM calculations to simulate the absorption (or the emission) spectra. The 100 snapshots have been extracted from the last 5 ns of the MD simulation (from the geometry stabilised region) by selecting one snapshot each 0.05 ns, leading to 100 snapshots. (2) In order to simulate these spectra up to 300 nm, seven electronic excited states have been considered for computing the excitation (or de-excitation) energies. In addition, the absorption (or emission) spectra have been simulated as a Gaussian convolution of the excitation (or de-excitation) energies of 100 snapshots, using a half-width at half-maximum (HWHM) of 0.2 eV.

Regarding the force field of the oxyluciferin forms (and their synthetic analogues), two sets of parameters have been obtained for each compound: one for the ground state (GS) and another for the excited state (ES). First, for each electronic state (GS and ES), a preliminary MD simulation has been performed with a set of parameters (bond distances, angles, dihedrals and Mulliken atomic charges) obtained in the gas phase for the ground minimum, and for the first excited state minimum respectively. One snapshot has been extracted from the geometry stabilised region (snapshot with the lowest energy) of the MD simulation and optimised at the QM/MM level of theory (at the GS, resp. at the ES). Second, a new set of parameters obtained from the snapshot optimised at the QM/MM level was extracted. A second MD simulation has been performed with this new set of parameters, respectively for the GS and ES. Again, one snapshot with low energy has been extracted from each MD simulation and optimised at the QM/MM level of theory. For each MD simulation with a particular set of force field parameters, a spectrum (absorption or fluorescent spectrum respectively) was simulated from 100 snapshots. From the 10 ns MD simulation, 100 snapshots have been extracted, and the QM/MM absorption (or emission) spectrum has been simulated (Figure 35). The iteration was stopped once two consecutives spectral shapes obtained with two consecutives parameters sets were converged. Usually, 2 or 3 iterations are enough for convergence.

The parameters achieving a good convergence between the simulated spectra have been applied for the final MD simulations.

(2) For more information regarding the extraction of the snapshots, see the section: Tutorial for System preparation -Input files given in the appendix part. Furthermore to check the spectral shape convergence of the simulated spectra with 100 snapshots, a spectrum of one compound has been simulated considering 1000 snapshots and the result shows that 100 snapshots were enough to reach the convergence (Figure 36a). Finally to check the suitability of the compound's parameters, absorption (emission) spectra simulated with and without optimising the compounds at the QM/MM level of theory have been performed. In all cases, similar spectral shapes have been obtained (Figure 36b). In this Ph.D. Thesis, the crystallographic structure of engineered firefly's luciferase (North American firefly, Photinus pyralis) has been downloaded from the Protein Data Bank (PDB code 4G37). 125 This firefly's luciferase is characterised by a disulfur bridge (linkage) between two residues, ILE108 and TYR447, which is built to keep the closed conformation of this luciferase, given the so-called 4G37 luciferase conformation. 125 The choice of using the 4G37 luciferase conformation (Figure 37) instead the open luciferase conformation (PDB code 4G36) 125 for this work was related to the fact that, during the process of bioluminescence, light emission occurs in this 4G37 luciferase conformation. 126

Classical molecular dynamics (MD) simulations: compound in protein

Figure 37 Presentation of firefly's luciferase in its closed conformation (PDB code 4G37).

The ligand (substrate) inside the crystallographic 4G37 luciferase conformation, 5'-O-[Ndehydroluciferyl-sulfamoyl]-adenosine (DLSA), has been replaced by the oxyluciferin structure (or by its analogues) and adenosine monophosphate (mono-deprotonated and further donated as AMPH) (Figure 38). In addition, the linkage between ILE108 and TYR447 residues has been deleted, as there is no accurate parameters to describe it, which are essential for performing the MD simulations. It was checked that the 4G37 luciferase conformation was held along the MD simulations although the linkage had been removed by computing the distance between the carbon alpha atom of the ILE108 and TYR447 residues along the MD simulation varying less than 2 Å.

P a g e | 83

Furthermore, the residues of the 4G37 luciferase conformation have been protonated by using LEaP from the Amber14 program. 123 In order to neutralise the system, seven histidine residues (76, 171, 310, 332, 419, 461 and 489) have been protonated after computing their pKa with the H++ program. 127 Each compound has been solvated using TIP3P water molecules model 124 within truncated octahedral box, ensuring a solvation shell of 10 Å around the compound, corresponding to around 19000 solvent molecules.

The force field parameters used for AMPH have been designed previously by our group. 4,[START_REF] Navizet | Color-Tuning Mechanism of Firefly Investigated by Multi Configurational Perturbation Method[END_REF]128 Amberff14SB force field parameters for have been used for the protein residues.

The force field parameters of the chemical compounds, their MD simulations of 20 ns and their absorption (or emission) spectra in protein environment have been performed

following the same procedure as in the case of the compounds in explicit water (see the section: compound in explicit water).

QM calculations

Density functional theory (DFT) has been used for the optimisation of the ground state (GS) minima. Time-dependent density functional theory (TD-DFT) was employed for the excited state (ES) optimisation and for the vertical transition energy calculations. The B3LYP 129,130 functional and the 6-311g(2d,p) basis set have been used for all calculations mentioned in the Chapters 6, 7 and 9. It have been previously reported that this functional and basis set are suitable for the compounds under study. [131][132][133][134] A benchmark with different DFT exchange-correlation functionals and basis sets has been previously achieved (see Tables S1-S3 reported in ref. 135 ). It should be noted that the B3LYP functional is reliable and suitable for the calculations as the compounds under study are characterised by small charge-transfer character (CT) of the transitions. While, when studying compounds with strong CT character, this functional usually yields un-physical results and other functional should be employed like the hybrid exchange-correlation functional.

The absorption (S0→S1 at the GS optimised geometry) and emission (S1→S0 at the first ES optimised geometry) energies have been computed in the gas phase and by considering implicit and explicit solvent (by using different solvation models). In all cases water molecules have been chosen as the solvent in line with the experiment. For implicit water solvent, the polarisable continuum model (PCM) 136 with a dielectric constant of 78.35 has been used. The compound structures have been optimised using the equilibrium regime within the linear P a g e | 84

response-PCM method. The vertical transitions have been performed using the non-equilibrium regime within the linear response (LR), corrected linear response (cLR) and state specific (SS) methods. The cLR and SS methods have been used only in the Chapter 7.

For checking the charge transfer (CT) character of the vertical transitions, the amount of charge transferred (q CT ) from one part to the other part of the molecule has been computed.

Finally, the natural transition orbitals (NTO) for the compounds have been plotted with Avogadro 137 package (see Tables S5 andS6 reported in ref. 135 ). All calculations mentioned above have been achieved by using the Gaussian 09 package. 138

QM/MM calculations

One hundred snapshots have been extracted along the classical MD simulation of each compounds and taken as the starting structures for simulating the absorption and emission spectra at the QM/MM level of theory with a non-polarisable scheme. All compounds (ligands)

have been described at the QM level while the water molecules, AMPH and the protein residues have been treated at the MM level within a non-polarisable model.

For each extracted snapshot, water molecules within a radius of 15 Å from the QM subsystem (ligand) have been selected and only the ligand and the surrounding protein residues within a radius of 9 Å were allowed to move. A coupling scheme 139 between Gaussian09 138 and Tinker 140 has been applied to perform the QM/MM calculations.

The QM/MM optimisations of the GS and the first ES have been performed, and then the S0→S1 and S1→S0 vertical transition energy have been computed, corresponding to the absorption and emission transitions, respectively. The functional B3LYP 129,130 and the basis set 6-311g(2d,p) have been used.

The electrostatic potential fitted (ESPF) 139 approach has been applied to compute the interaction between the charge density of QM subsystem and the external electrostatic potential of the MM subsystem within 9 Å from the QM subsystem. ESPF approach permits to take into account the external electrostatic potential created by the MM point charges in the calculation of the QM wavefunction. The micro-iterations technique 141 has been used to converge the MM subsystem for every minimisation step of the QM subsystem.

For the absorption and emission spectra simulation, a Gaussian convolution of the computed vertical transition energies (corresponding to the absorption and emission P a g e | 85

transitions, respectively) of 100 snapshots has been performed, using a half-width at halfmaximum (HWHM) of 0.2 eV.

Vibronic spectra

On the basis of the TD-DFT vibrational signatures, the vibronic spectra have been generated by using the Franck-Condon (FC) and the Franck-Condon-Herzberg-Teller (FCHT)

approximations. These FC and FCHT approximations are implemented in Gaussian 09. 109,111,121,122 The absorption and emission vibronic spectra have been generated in gas phase and in solution considering different solvation models (PCM model and microsolvation coupled with PCM model). For the microsolvation coupled with PCM model, 5 snapshots have been extracted

from the MD simulation and only water molecules within a radius of 3 Å around the heteroatoms of the ligand have been kept. Then, these snapshots have been optimised in the GS using the same level of theory, functional and basis set as mentioned above (see QM calculations).

Afterwards, they have been optimised in the ES keeping frozen the water molecules at the same orientation of the GS minima, because the geometry of the GS and ES minima should be similar when generating the vibronic spectra within a harmonic adiabatic approximation.

The band broadening of the vibronic spectra has been simulated by a Gaussian function with HWHM of 135 cm -1 . For a comparative purpose, the absorption and the emission vibronic spectra have been shifted to match the spectra maxima obtained with the PCM model, that is, they have been shifted to 2.17 eV for the absorption and to 2.16 eV for emission.

Regarding calculations mentioned in Chapter 8, the same procedure as detailed above has been applied to compute the vibronic absorption spectra and the vertical transition energies. For simulating the vibronic spectra we applied the Vertical Gradient (VG) model instead the Adiabatic Hessian (AH) one as the relative positions of the adenosine monophosphate (AMP) and oxyluciferin analogue are not the same in the ground and excited states (AH calculations are infeasible due to the flexible complexes including the analogue form and AMP). All calculations have been achieved by using the Gaussian 16 package. 142 P a g e | 86

Chapter 6

Oxyluciferin and its synthetic analogues in water solution

Context

During the firefly's bioluminescence process, the molecule responsible for light emission, oxyluciferin, can take several chemical forms due to inter-exchange reactions such as keto/enol tautomerization and the phenol and/or enol deprotonation/protonation equilibria (Figure 2). This triple equilibrium makes the oxyluciferin forms sensitive to the solvent, pH changes, specific interaction with bases and polarity of the luciferase active site. All these factors complicate the understanding of the bioluminescence process, and in particular the oxyluciferin chemical form contributing to the light emission. A theoretical investigation study on the oxyluciferin chemical forms in different solvents and considering the protein environment has been published in 2011 by Chen et al. [START_REF] Chen | Systematic theoretical investigation on the light emitter of firefly[END_REF] The results of this study exclude all the neutral oxyluciferin forms as light emitters while, anionic oxyluciferin forms are considered as emissive forms and in particular, the phenolate-keto form. Despite the numerous theoretical and experimental studies completed throughout the last decade, [START_REF] Navizet | Color-Tuning Mechanism of Firefly Investigated by Multi Configurational Perturbation Method[END_REF][16][17][START_REF] Ugarova | Protein structure and bioluminescent spectra for firefly bioluminescence[END_REF][START_REF] Hirano | Spectroscopic Studies of the Light-Color Modulation Mechanism of Firefly (Beetle) Bioluminescence[END_REF] the debate on the chemical nature of the light emitter is still ongoing.

The purpose of this Chapter is to demonstrate whether oxyluciferin synthetic analogues (Figure 3) can be used to mimic the natural oxyluciferin forms in water solution. The analogue forms avoid the inter-exchange reactions, that is, the triple equilibria (Figure 2). The oxyluciferin forms and their analogues studied in this Chapter are reported in Figure 3 which contains 12 compounds: 6 forms of oxyluciferin and 6 forms of their analogues, including the protonated (neutral) and deprotonated (anionic) forms. The phenol-cycle (5,5-spirocyclopropyloxyluciferin) analogue which has a cyclopropyl group in the thiazolone moiety is designed to block the ketoenol tautomerization equilibrium. The phenol-OMe (4-methoxyoxyluciferin) analogue which contains a methoxy group in the thiazolone moiety blocks not only the keto-enol tautomerization, but also the enol deprotonation equilibrium. The OMe-enol (6′methoxyoxyluciferin) analogue which has a methoxy group in the benzothiazole moiety blocks the phenol deprotonation equilibrium.

P a g e | 87

For this aim, our collaborators have measured experimentally the absorption and emission spectra of three oxyluciferin analogues (phenol-cycle, phenol-OMe and OMe-enol forms in Figure 3) in water at pH 5 and 11 using the steady-state spectroscopy experiments. In parallel to this experimental measurements, we have simulated the absorption and emission spectra of each oxyluciferin analogue (protonated and deprotonated forms) and its corresponding natural form (Figure 3) in explicit water molecules, taking into account the solvent-solute interactions and dynamical effect of the system by using classical MD simulations and hybrid QM/MM methods. To realise this work, we have performed a classical MD simulation for each compound in explicit water molecules in order to obtain a statistical number of snapshots. Afterward, we have carried out the QM/MM calculations to compute the vertical transition energies and then simulated the absorption and the emission spectra (see computational details in Chapter 5). We have compared the measured and simulated spectra of the analogues with the ones of their corresponding oxyluciferin natural forms to give insight into the analogues ability to mimic oxyluciferin forms in water.

A quite good agreement between the simulated and experimental spectra of the oxyluciferin analogues in water has been achieved. From this comparative study of the photochemical properties of the oxyluciferin natural forms and the synthetic analogues, we demonstrate that the oxyluciferin analogues are good candidates for mimicking natural oxyluciferin forms in water solution.

Results and discussion

H-Bond interactions between Solute-Solvent

Before studying the photochemical properties of each compound, we have analysed the H-bonding formed between each compound (the oxyluciferin and analogue forms) and water molecules during the classical MD simulations. Indeed, it has been reported in previous studies that these interactions impact both the absorption and emission spectra. 131,143,144 For this, we have calculated using the CPPTRAJ package the statistical average number of H-bonding between the heteroatoms (O10, O11, N3 and N3′) of each compound and water molecules during the MD simulations in the excited state (ES) (Figure 39D). 39A, B and C, respectively), the deprotonated oxygen (O10 and O11) shows 3 to 4 H-bond interactions involving 3 or 4 different water molecules. On the other hand, when this oxygen atom of the compound structure is protonated, it interacts with only one water molecule forming one H-bond interaction (oxygen of the enol group in Figure 39B or phenol group in Figure 39C). Moreover, each nitrogen atom of the compound structures shows in general one H-bond interaction during the MD simulation.

Furthermore, similar average number of H-bonds was observed for the oxyluciferin analogues (Figure 39D) indicating that the number of H-bond interactions is impacted neither by the change to cyclopropyl nor the change to methyl substituents. Only a slight difference was found for phenolate-OMe and OMe-enolate compared to their corresponding oxyluciferin forms, phenolate-enol and phenol-enolate, respectively, which can be caused by the inductive effect of the methyl group that increases the charge of the oxygen atom. The fact that the introduction of the substituents (cyclopropyl and methyl groups) within the oxyluciferin structure does not significantly modify the H-bond networks along the MD simulation is a good sign indicating that the analogue forms are able to imitate their corresponding oxyluciferin forms in water.

Figure 39 Graphical illustration of the H-bond pattern of representative MD snapshots of A) phenolateketo, B) phenolate-enol, and C) phenol-enolate (Left side). D) Average number of H-bond interactions between the heteroatoms (O10, O11, N3 and N3′) of the compound and water molecules found along the MD trajectory (Right side). Each oxyluciferin form is depicted in the same colour as its corresponding

analogue (Figure reproduced from ref. [START_REF] Garcia-Iriepa | Simulation and Analysis of the Spectroscopic Properties of Oxyluciferin and Its Analogues in Water[END_REF] ).

P a g e | 89

Simulation of the absorption spectra

We have compared the experimental absorption spectra obtained with the synthetic analogues (measured by our collaborators) with the ones simulated using QM/MM calculations.

The simulated spectra have been obtained as a Gaussian convolution of excitation energies of 100 MD snapshots computed at the QM/MM level of theory considering 7 electronic excited states. For comparative purpose the simulated absorption spectra are given as a normalised intensity against wavelengths in nanometer (nm) to match the ones measured experimentally.

The simulated and the experimental absorption spectra of the analogues are reported in Figure 40. Apart from a small difference found in term of broad spectral shape (experimental spectra are broadest than the simulated ones), a good agreement between the simulated and experimental spectra was obtained in terms of both the wavelength absorption maxima and the spectral shape (Figure 40). The narrowing of the simulated spectra may be due to the nonconsideration of the vibronic coupling in the calculations, which in some cases [START_REF] Jacquemin | Impact of Vibronic Couplings on Perceived Colors: Two Anthraquinones as a Working Example[END_REF][START_REF] Charaf-Eddin | Vibronic Spectra of Organic Electronic Chromophores[END_REF] has a significant influence on the spectra simulation (This hypothesis will be discussed in the next Chapter). Furthermore, in order to provide a general overview of the ability of different analogues to mimic their corresponding natural forms, we compared their simulated absorption spectra in water (Figure 41). In general, a good agreement between the simulated absorption spectra of the analogues and the ones of their corresponding natural forms has been achieved. Except for the phenol-cycle analogue, which presents a dominating absorption band slightly different from the one simulated for the phenol-keto (Figure 41C). Moreover, we compared the simulated absorption spectra of the oxyluciferin forms in water with the reported ones measured by applying the multivariate curve resolution-alternating least-squares (MCR-ALS) procedure. 16 P a g e | 90

Apart from the small energy differences between the simulated and experimental spectrum maxima (less than 0.1 eV), all the simulated spectra are in line with the experimental ones.

Analysis of the involved electronic transition

Additionally to the spectral shapes of the absorption spectra, we have analysed the electronic nature of the absorption transitions of the compound structures. For instance in the case of phenol-keto and its analogue, which exhibit one of the broadest absorption bands (Figure 41C), the computed first three singlet (S1, S2, and S3) excited states are optically bright and close in energy so that the dominant absorption band corresponds not only to the electronic transition S0→S1 with the low-energy, but also to S0→S2, and S0→S3 electronic transitions, explaining the broad and structured spectral shape of the simulated and experimental spectra (Figure 40A and Figure 41C). In order to analyse the electronic nature of these three electronic transitions, we investigated the molecular orbitals involved in each transition (Figure 42). For the phenol-keto and its analogue, the electron density delocalises from the benzothiazole moiety to the thiazolone one during the S0→S1 transition. The S0→S2 transition is similar to the S0→S1 one, but some of the electron density of the bonding orbital, HOMO-1, are located on P a g e | 91

the sulphur (S) atom of the benzothiazole moiety (Figure 42). In addition, S1 and S2 states are very close in energy, and their electronic nature results in a mixture of HOMO-LUMO and HOMO-1-LUMO mono-excitations, which can explain the similarity of their absorption intensities.

Regarding the phenol-enol and its two analogues (deprotonated and deprotonated forms), the S1 and S2 states are energetically more distant than the ones obtained with phenolketo (or with its analogue) leading to two separated absorption bands (Figure 41). Moreover, the electronic nature of S0→S1 and S0→S2 transitions is similar to the ones observed for the phenol-keto and its analogue form (Figure 42 and . In all cases, the S0→S1 transition corresponds to the band absorption maximum which is more intense than the S0→S2 transition, in particular the S0→S2 transitions of phenolenol and its analogues (Figure 41A andE).

Hence, we can conclude that the electronic nature of the vertical transitions and absorption spectral band of the oxyluciferin analogues are almost identical to that of the natural oxyluciferin forms in water, showing that these analogues are good candidates for mimicking the oxyluciferin natural forms in water.

Simulation of the emission spectra

The experimental emission spectra of three oxyluciferin analogues (phenol-cycle, phenol-OMe and OMe-enol forms in Figure 3) in water at pH 5 and 11 have been measured by our collaborators using the steady-state spectroscopy. The experimental group has observed that the emission spectra recorded at pH 5 perfectly overlap the ones obtained at pH 11, except for the phenol-OMe analogue which shows at pH 5 two distinct bands, an intense band around P a g e | 92

550 nm and weaker one around 450 nm (dashed blue and red lines Figure 43). In parallel to this experiment, we have simulated the emission spectra of each oxyluciferin analogue (protonated and deprotonated forms) and its corresponding natural form (Figure 3) using QM/MM methods considering explicitly the water molecules. By comparing the recorded spectra of the analogues at pH 5 and 11 with the simulated ones of the protonated and deprotonated analogue forms, we found that the simulated spectra of the deprotonated analogues are in good agreement with experimental ones (at pH 11) in terms of both the wavelength emission maxima and the spectral shape (dashed and continuous blue lines in Figure 43). While, the simulated spectra of the protonated analogue forms do not overlap those of the deprotonated forms, contrary to what has been observed experimentally (continuous blue and red lines in Figure 43). The disagreement between the experimental and the simulated spectra of protonated analogues forms is due to the very fast excited state proton transfer (ESPT), which prevents the detection of the emission from the protonated forms, and hence only the emission of the deprotonated analogues (at pH 11) is recorded by the steady-state fluorescence experiment.

In order to demonstrate the ultrafast ESPT, the experimental group has performed femtosecond pump-probe measurements on phenol-cycle at pH 5 (protonated) and at pH 11 Furthermore, we have also simulated the emission spectrum for each natural oxyluciferin form and compared it with the one simulated for the corresponding analogue (Figure 44). In general, the energy differences between the simulated emission maxima of each P a g e | 93

oxyluciferin form and its analogue are less than 0.12 eV. The emission spectra of phenolate-enol and phenolate-keto are slightly blue-shifted with respect to the spectra of their analogues, while the emission spectrum of phenol-enol is slightly redshifted (Figure 44). The emission spectra of phenol-enolate and phenol-keto match the spectra of their corresponding analogues.

Moreover, we have compared the simulated emission spectra of the oxyluciferin forms in water with the reported ones measured by applying the MCR-ALS procedure. 16 Apart from the small energy differences between the simulated and experimental spectrum maxima (less than 0.1 eV), all the simulated spectra are in line with the experimental ones.

Conclusion

In this work, we have simulated the absorption/emission spectra of different oxyluciferin forms and their analogues in explicit water taking into account the solvent-solute interactions and dynamical effect of the system. In addition, we have compared the simulated spectra with the ones measured by the experimental group.

The conclusion to be drawn from this work is that the electronic structures and absorption/emission spectra of the oxyluciferin analogues under study are quite similar to their P a g e | 94

corresponding natural forms in water solution. The introduction of the cyclopropyl and methyl groups within the oxyluciferin structure does not change the spectral shape nor the wavelength maxima of the absorption/emission spectra.

In general, a quite good agreement between the simulated and steady-state fluorescence spectra of the oxyluciferin analogues in water has been achieved. In addition, as the spectra of the protonated compounds cannot be recorded by using steady-state experiment, the simulated spectra of oxyluciferin natural forms in water have been compared and are consistent with the reported ones measured by applying the MCR-ALS procedure. 16 Moreover, the experimental group has demonstrated that, by using femtosecond pump-probe measurements on phenol-cycle at pH 5 and 11, the protonated phenol-cycle form is ruled by an ultrafast ESPT reaction, which can explain the non-recording of its spectrum with a steady-state experiment.

Finally from the comparative study of the photochemical properties of all the compounds under study, we can state that the oxyluciferin synthetic analogues are able to mimic natural oxyluciferin forms in water solution. The calculation of vertical transition energies is a primordial step to simulate the absorption and emission spectra of a given molecular system. However, the calculation of these transition energies is not always an easy task, especially when the surrounding effects (e.g., solvent or protein environment) have to be taken into account. The difficulty of these calculations comes from the interaction between the relaxation processes of both the solute and the solvent. There are various solvation models that describe differently the solvent effect on the electronic structure of the solute. For instance, the explicit solvation model maintains an atomistic representation of the solvent molecules and hence, it considers the interactions between solute-solvent molecules allowing a more accurate simulation of the environment.

However, this model requires a large number of solvent conformations to achieve optimal vertical electronic transition energies and simulate complete absorption and emission spectra.

One the other hand, other solvation models such as the implicit PCM model [START_REF] Tomasi | Quantum mechanical continuum solvation models[END_REF] describe the solvent as a polarisable dielectric characterised by a dielectric constant (ε) and therefore it does not take into account the explicit solute-solvent interactions. In the framework of excited-state (ES) calculations within the PCM model, the response of solvent to the change of the solute charge density is taken into account differently depending on the chosen formalism: linear response (LR), [START_REF] Cossi | Time-dependent density functional theory for molecules in liquid solutions[END_REF][START_REF] Cammi | Fast evaluation of geometries and properties of P a g e | 205 excited molecules in solution: a Tamm-Dancoff model with application to 4dimethylaminobenzonitrile[END_REF] corrected linear response (cLR) [START_REF] Caricato | Formation and relaxation of excited states in solution: A new time dependent polarizable continuum model based on time dependent density functional theory[END_REF] and state specific (SS) [START_REF] Cossi | Solvent effect on vertical electronic transitions by the polarizable continuum model[END_REF] formalisms.

The purpose of this Chapter is to highlight and evaluate the solvation effect on the computed vertical transition energies and hence on the simulated absorption and emission spectra of oxyluciferin natural forms and the synthesised analogues (Figure 3). For this aim, we used both the explicit and the implicit solvation models. For the calculation with explicit model, we performed a classical MD simulation for each compound in explicit water molecules in order to obtain a statistical number of snapshots. Afterward, we carried out the QM/MM calculations to compute the vertical transition energies and then simulate the absorption and the emission P a g e | 106

spectra. For the calculation with implicit model (PCM), we computed the transition energies within three formalisms, LR, cLR and SS formalisms, at the QM level of theory (see computational details in Chapter 5). Then, we analysed the spectral shape of the simulated spectra and compared their maximum wavelengths to the ones measured experimentally. 16,[START_REF] Garcia-Iriepa | Simulation and Analysis of the Spectroscopic Properties of Oxyluciferin and Its Analogues in Water[END_REF] Finally, as considering the vibronic information in the calculations could be essential to reproduce the overall spectral shape, [START_REF] Jacquemin | Impact of Vibronic Couplings on Perceived Colors: Two Anthraquinones as a Working Example[END_REF][START_REF] Charaf-Eddin | Vibronic Spectra of Organic Electronic Chromophores[END_REF][START_REF] Santoro | Going beyond the vertical approximation with timedependent density functional theory[END_REF][START_REF] Ferrer | Insights for an accurate comparison of computational data to experimental absorption and emission spectra: beyond the vertical transition approximation[END_REF][START_REF] Pedone | Environmental effects in computational spectroscopy: accuracy and interpretation[END_REF][START_REF] Bloino | Aiming at an accurate prediction of vibrational and electronic spectra for mediumto-large molecules: an overview[END_REF] we simulated the vibronic spectra of phenolate-keto oxyluciferin form in Gas Phase (GP), in PCM and with a microsolvation model (considering explicitly some solvent molecules around the solute in both GP and PCM) using the Franck-Condon (FC) 121,145 and Franck-Condon-Herzberg-Teller (FCHT) 111 approximations. (3) In this study, we provide a general analysis of the solvation model effects on the simulation of the absorption and emission spectra. We demonstrate that, the implicit SS-PCM model is the most suitable to compute the vertical transition energies of the oxyluciferin forms with a large charge transfer (CT) character while, the implicit LR-PCM model could be used for the chemical forms characterised by a small (almost negligible) CT character. Moreover, we demonstrated that the use of an explicit solvation model is needed to reproduce the measured experimental spectra shape. Finally, the vibronic spectra using an implicit solvation model show some noticeable differences compared to the ones in gas phase, while small changes were found when adding some explicit water molecules around the solute within a microsolvated model.

Results and discussion

Vertical transition energies

The S0→S1 and S1→S0 vertical transition energies, corresponding to the absorption and emission transitions respectively, of the oxyluciferin natural forms and their respective analogues are computed in gas phase and in both implicit (using LR-PCM, cLR-PCM and SS-PCM formalisms) and explicit solvation models. The computed vertical transition energies, the corresponding wavelength spectra maxima and the experimental results for the absorption and emission transition are regrouped in Table 1 and Table 2, respectively.

(3) FC approximation considers that the electronic transition dipole moment is constant during the electronic transition (independent of the nuclear coordinates around the initial and final equilibrium geometries). HT approximation considers a small variation of the relative positions and velocities of the nuclei during the electronic transition.

Absorption transition energies

If we compare the absorption transition energies computed using the LR, cLR and SS formalisms, we observe that the absorption energies are mostly underestimates compared to experimental ones, except for the energies of phenolate-keto form and its analogue which are slightly above the experimental ones, being the energy difference of less than 0.06 eV. The energy differences calculated between LR-PCM (𝐸 𝐿𝑅 𝑎𝑏𝑠 ) and cLR-PCM (𝐸 𝑐𝐿𝑅 𝑎𝑏𝑠 ) are small (≤ 0.1 eV).

However, the transition energies computed using SS-PCM (𝐸 𝑆𝑆 𝑎𝑏𝑠 ) are significantly different from the ones computed with LR-PCM and cLR-PCM, except for phenolate-keto, phenol-enol and their analogues.

In order to rationalise the energy differences found between LR-PCM (𝐸 𝐿𝑅 𝑎𝑏𝑠 ) and SS-PCM (𝐸 𝑆𝑆 𝑎𝑏𝑠 ), we computed the charge transfer (CT) character of the absorption transition (S0→S1) for the oxyluciferin natural forms and the analogues as it is known that the LR formalism is suitable for transitions characterised by low CT character, unlike the SS formalism which gives better results for transitions with large charge rearrangement. 146,147 In particular, we calculated the Mulliken charges of the thiazolone moiety (qS0 thiaz and qS1 thiaz ) for the GS equilibrium geometry at the initial (S0) and final (S1) states, respectively. Then, we calculated the CT character as the difference of the Mulliken charges of the thiazolone moiety between the final state qS1 thiaz and initial state qS0 thiaz (CT = qS1 thiaz -qS0 thiaz ) (Table 1). By analysing the calculated CT character for the absorption transition ( Regarding LCT group, we observed that absorption energies computed within LR, cLR and SS formalisms are quite similar presenting energy differences less than 0.1 eV (Table 1). This can be explained by the negligible CT character of the absorption transition. For the HCT group, the energy differences between 𝐸 𝐿𝑅 𝑎𝑏𝑠 and 𝐸 𝑆𝑆 𝑎𝑏𝑠 calculated for the phenol-keto and phenolenolate forms and their corresponding analogues are the largest ones, about 0.27-0.37 eV in P a g e | 108

absolute value (Table 1) due to the large CT character of the transition. The absorption energy of the phenol-keto and phenol-enolate forms and their corresponding analogues (HCT group) is more sensitive to the solvation model. For the MCT group, we observe a halfway effect of the formalism chosen, finding out that the transition energies computed with SS formalism are always higher than the ones computed with LR formalism (energy difference around 0.13 to 0.24 eV).

Table 1 Absorption energies (4) (given in eV (nm)) of oxyluciferin natural forms (coloured lines) and their analogues (uncoloured lines) computed in gas phase (GP), using implicit (LR-, cLR-and SS-PCM) and explicit (5) (QM/MM) solvation models. Experimental absorption wavelength maxima (given in eV (nm))

from ref. 16 are also presented. CT character is given as CT = qS1 thiaz -qS0 thiaz . Minus sign indicates that, during the transition, there is a fraction of the electronic charge transferred from the benzothiazole moiety to the thiazolone moiety. (Table reproduced from ref. 135 ).

Furthermore, we compared the absorption energies computed within LR, cLR and SS formalisms to the experimental results (Table 1). As previously mentioned for phenol-enol and its analogue (LCT group), the absorption energies computed with the LR, cLR and SS formalisms are quite similar. For LCT group, the absorption energies (𝐸 𝑐𝐿𝑅 𝑎𝑏𝑠 ) computed using cLR formalism is the closest to the experimental (𝐸 𝜆𝑚𝑎𝑥,𝑒𝑥𝑝

𝑎𝑏𝑠

) values, with an energy difference about 0.05 eV.

Then for HCT and MCT groups where the CT character of the chemical compounds are large, the (4) The transition energies have been computed using the B3LYP functional and the 6-311g(2d,p) basis set.

A benchmark with different DFT exchange-correlation functionals has been achieved. The obtained results show that the B3LYP functional gives transition energies closest to the experimental values within a shorter CPU time (Table A-1). (5) 𝐸 𝑄𝑀𝑀𝑀 𝑎𝑏𝑠 values correspond to the maxima wavelengths of the absorption spectra simulated considering 100 MD snapshots.

P a g e | 109

absorption energies (𝐸 𝑐𝐿𝑅 𝑎𝑏𝑠 and 𝐸 𝑆𝑆 𝑎𝑏𝑠 ) computed using both cLR and SS formalisms are closest to the experimental values (Table 1). Except for phenol-keto and its analogue which constitute a special case. In this case, the closest value to the experiment is obtained by when using LR formalism with an energy difference less than 0.08 eV. This can be explained by the fact that the simulated and experimental absorption spectra of phenol-keto (or of its analogue) are different from the ones obtained for the rest of chemical compounds. For most of chemical compounds, the absorption spectra are characterised by an intense band that corresponds to the S1 transition (described by a HOMO to LUMO mono-excitation) and by a less intense one that corresponds to S2 transition. However, the absorption spectral band of phenol-keto and its analogue is wide and structured because both S1 and S2 transitions are intense and close in energy, as we have discussed in the previous work (Chapter 6). In addition, the electronic nature of these S1 and S2

transitions results in a mixture of mono-excitations: HOMO to LUMO and HOMO-1 to LUMO.

Moreover, we also compared the absorption energies computed using explicit solvation model to the experimental results. For all the oxyluciferin forms and their analogues, the absorption energies [START_REF] Mochizuki | Robust red-emission spectra and yields in firefly bioluminescence against temperature changes[END_REF] (𝐸 𝑄𝑀𝑀𝑀 𝑎𝑏𝑠

) computed using explicit solvation model match the experimental values, with an energy difference less than 0.1 eV (Table 1).

In summary, although different formalisms have been applied to account the solvent effect within PCM solvation model, the absorption transition energies of some of the oxyluciferin forms and their analogues are still far from the experimental values. In addition, the computed transition energies differ from one formalism to another, depending on the CT character of the chemical compound. For instance, the SS-PCM formalism is suitable to compute the absorption energies of the chemical compounds characterised by a significant CT character while, the LR-PCM formalism is appropriate for the compounds with negligible CT character.

Furthermore, the lack of explicit interactions between the solute and solvent molecules could also be the reason for the discrepancy between these energies. Indeed, the vertical transition energies computed considering explicit solvent model are in good agreement with experimental values. [START_REF] Mochizuki | Robust red-emission spectra and yields in firefly bioluminescence against temperature changes[END_REF] It is worthwhile to highlight that in order to explore a high enough conformational space it would have been desirable to perform thousands of conformations (snapshots) to simulate the absorption and emission spectra. This is due to the fact when considering the explicit water molecules in the system the degrees of freedom increase very rapidly. Because of the physical constraints (such as time and computational resources), we could only use one thousand snapshots. Looking further into the spectra convergence, we were able to demonstrate that with only one hundred snapshots reasonably accurate spectra can be simulated. Indeed, the spectra simulated with one thousand and with one hundred snapshots are similar (see Figure 36a in computational details).

Emission transition energies

In the following, we performed a similar reasoning for the emission transition energies.

If we compare the emission transition energies computed using the LR, cLR and SS formalisms (Table 2), we observe that the emission energies (𝐸 𝐿𝑅 𝑒𝑚𝑖 and 𝐸 𝑐𝐿𝑅 𝑒𝑚𝑖 ) computed with LR and cLR formalisms are similar, showing energy differences less than 0.1 eV as found for the absorption transition, except for phenol-keto and its analogue. We also analysed the CT character of the emission transition (S1→S0) and we calculated the charge difference of the thiazolone moiety at the ES equilibrium geometry between the S1 and S0 states (CT = qS0 thiaz -qS1 thiaz ) (Table 2). The CT = qS0 thiaz -qS1 thiaz is similar to its opposite CT = qS1 thiaz -qS0 thia obtained for the absorption transition and hence, we kept the same LCT, MCT and HCT groups to discuss the CT character of the emission transition.

Table 2 Emission energies (given in eV (nm)) of oxyluciferin natural forms (coloured lines) and their

analogues (uncoloured lines) computed in gas phase [START_REF] Mofford | Aminoluciferins Extend Firefly Luciferase Bioluminescence into the Near-Infrared and Can Be Preferred Substrates over D-Luciferin[END_REF] (GP), using implicit (LR-, cLR-and SS-PCM) and explicit [START_REF] Reddy | Robust Light Emission from Cyclic Alkylaminoluciferin Substrates for Firefly Luciferase[END_REF] (QM/MM) solvation models. Experimental absorption wavelength maxima (given in eV (nm)) from ref. 16 are also presented. CT character is given as CT = qS0 thiaz -qS1 thiaz . Minus sign indicates that, during the transition, there is a fraction of the electronic charge transferred from the thiazolone moiety to the benzothiazole moiety. (Table reproduced from ref. 135 ). [START_REF] Mofford | Aminoluciferins Extend Firefly Luciferase Bioluminescence into the Near-Infrared and Can Be Preferred Substrates over D-Luciferin[END_REF] "For phenol-keto and its analogue in GP, the electronic nature of S1, after optimisation in this state, is not the expected π,π* state but a n,π* dark state with the n orbital centred on the oxygen atom of the keto moiety. In fact, the n,π* state is over-stabilised in GP. For this reason, we do not give their emission energy in GP. However, the right ordering of the ES is recovered when considering the solvent effect, either implicitly or explicitly". 135 [START_REF] Reddy | Robust Light Emission from Cyclic Alkylaminoluciferin Substrates for Firefly Luciferase[END_REF] 𝐸 𝑄𝑀𝑀𝑀 𝑒𝑚𝑖 values correspond to the maxima wavelengths of the emission spectra simulated considering 100 MD snapshots.

P a g e | 111

For the LCT group, that is, the phenol-enol and its analogue, we observed that the emission energies computed within LR, cLR and SS formalisms are similar, presenting energy differences less than 0.1 eV as found for the absorption transition (Table 2). This can be explained by the negligible CT character of the emission transition. Regarding the HCT group (phenol-keto, phenol-enolate and their analogues) where the CT character is important, we found the largest energy differences between LR and SS formalisms. It should be noted that for the chemical compounds of the HCT group, the emission energies computed using SS formalism (𝐸 𝑆𝑆 𝑒𝑚𝑖 ) is the closest to the experimental values with an energy difference less than 0.1 eV (Table 2). For MCT group (phenolate-keto, phenolate-enol, phenolate-enolate and their analogues)

characterising by an intermediate CT character, we observe different behaviours without a significant emission energy change between LR, cLR and SS formalisms (less than 0.1 eV).

Furthermore the computation of the vertical transition energies, we calculated the Stokes shifts (𝐸 𝑎𝑏𝑠 -𝐸 𝑒𝑚𝑖 ) with the different solvation models (Table 3). The Stokes shifts calculated using the LR, cLR and SS formalisms are mostly underestimated comparing to experimental ones. This finding results from the general underestimation of the absorption transition energies and an overestimation of the emission ones, as already discussed above (Table 1 and Table 2). From the Table 3, we observe similar trends for the Stokes shifts as those observed for the vertical transition energies, that is: (i) in the case of LCT group, the Stokes shift and transition energy values computed using the LR, cLR and SS formalisms are similar and close to the experimental values (due to the negligible CT character), (ii) in the case of the HCT and MCT groups, the Stoke shift computed with the SS formalism is the closest to the experiment one but (iii) the differences found between the Stokes shifts computed using SS formalism and those computed with the LR and cLR are smaller for the compounds characterised by an intermediate CT character (MCT group), in comparison to the compounds with a larger CT character (HCT group). Moreover, we also compared the emission energies computed using explicit solvation model to the experimental results. For all the oxyluciferin forms and their analogues, the absorption energies (𝐸 𝑄𝑀𝑀𝑀 𝑒𝑚𝑖 ) computed using explicit solvation model match the experimental values, with an energy difference less than 0.1 eV (Table 2).

P a g e | 112

In summary, although different formalisms have been applied to account the solvent effect within PCM solvation model, the emission transition energies of some of the oxyluciferin forms and their analogues are still far from the experimental values. In addition, the computed transition energies differ from one formalism to another, depending on the CT character of the chemical compound. For instance, the SS-PCM formalism is suitable to compute the emission energies of the chemical compounds characterised by a significant CT character while, the LR-PCM formalism is appropriate for the compounds with negligible CT character. Furthermore, the lack of explicit interactions between the solute and solvent molecules could also be the reason for the discrepancy between these energies. Indeed, the vertical transition energies and the corresponding Stokes shifts computed considering explicit solvent model are in good agreement with experimental values. Therefore, in general the consideration of explicit solvent molecules and hence, considering their interactions with the solute, is primordial to reproduce the experimental energy values.

P a g e | 113

Absorption and emission spectral shape

Apart from the computation of the vertical transition energies, we also analysed the general shape of the absorption and emission spectra. In the following, we will study the effect of considering a statistical number of conformations (snapshots) and vibrational normal modes to simulate the absorption and emission spectra. For this aim, we compared the QM/MM simulated spectra (considering 100 snapshots) to those simulated in gas phase and with the LR-PCM formalism (Figure 45). In general, the spectra obtained from statistical QM/MM calculations are more structured than the ones from single point gas phase or PCM calculations, in particular the absorption spectra of phenol-keto (Figure 45a). In addition, the absorption bands computed considering 100 snapshots are wider than the computed ones in gas phase or in PCM. A consequence of this finding is that when a second optically bright excited state (Sx) is relatively close in energy to the first optically bright excited state (S1), such as in the case of phenol-enol, its corresponding band can be hidden in the tail of the most intense (S1) band (Figure 45b). This broadening of the absorption spectra is in line with the measured experimental spectra. 16 Moreover, for some oxyluciferin forms, we observe a clear change of the relative intensities of the different absorption bands simulated considering multiple number of snapshots. For instance, the absorption spectra of phenol-enolate simulated in both gas phase and in PCM do not show any energy band corresponding to the second optically bright excited state, as its intensity is quite low (Figure 45c). However, when simulating these spectra with the explicit model considering 100 snapshots, a second band appears with marked intensity, in line with the measured experimental spectra. 16 P a g e | 114 Regarding the simulated emission spectra, the main difference observed when considering 100 snapshots is that the emission spectral band is broader and not symmetric (Figure 46).

P a g e | 115

In order to study the effect of the solvent model nature considering the vibrational information in the spectra simulation, we simulated the vibronic absorption and emission spectra of phenolate-keto [START_REF] Wang | Impact of site-directed mutant luciferase on quantitative Green and Orange/Red Emission Intensities in Firefly Bioluminescence[END_REF] in gas phase, in LR-PCM model and using a microsolvation model (considering a small sphere of explicit water molecules around the solute) to include the H-bond interactions (Figure 47). To generate the absorption and emission vibronic spectra, we tested two approximations: Franck-Condon (FC) and Franck-Condon-Herzberg-Teller (FCHT)

approximations. The vibronic spectra of phenolate-keto computed with the FC approximation resemble to the ones generated using FCHT approximation, showing only small differences in the band intensity (Figure 47). In the following, we present separately the vibronic absorption and emission spectra and analyse the involved vibronic transitions computed with FCHT approximation (Figure 48a and 51b). From the vibronic absorption (Figure 48a) and emission (Figure 48b) spectra similar trends can be observed. First, the spectra computed in gas phase (GP) regardless if the microsolvation (µSolv) effect has been considered or not are broader than the ones obtained with the implicit model (PCM) model (Figure 48a and Figure 48b). In addition, we observe that the band following the most dominant band for the spectra computed in GP is more intense than those computed with PCM model (beige shadow in Figure 48 and Figure 48b). We note that the dominant band (grey shadow in Figure 48a and Figure 48b) corresponds to the 0-0' and 0'-0 vibronic transitions in all cases except for the one in not microsolvated GP. [START_REF] Wang | Impact of site-directed mutant luciferase on quantitative Green and Orange/Red Emission Intensities in Firefly Bioluminescence[END_REF] The choice of using the phenolate-keto form of oxyluciferin is based on assumption that it is the most probable chemical form responsible of the emission in fireflies bioluminescence. [START_REF] Chen | Systematic theoretical investigation on the light emitter of firefly[END_REF][START_REF] Naumov | Spectral-Structural Effects of the Keto-Enol-Enolate and Phenol Phenolate Equilibria of Oxyluciferin[END_REF][START_REF] Branchini | Yellow-green and red firefly bioluminescence from 5,5-dimethyloxyluciferin[END_REF][START_REF] Rebarz | Deciphering the protonation and tautomeric equilibria of firefly oxyluciferin by molecular engineering and multivariate curve resolution[END_REF][START_REF] Naumov | Spectral-Structural Effects of the Keto-Enol-Enolate and Phenol Phenolate Equilibria of Oxyluciferin[END_REF] P a g e | 116

Moreover, we can analyse the implicit solvation and microsolvation effects not only on the spectral shape, but also on the vibronic transitions. If we compare the not microsolvated systems in GP and with PCM model, we observe that all the medium-intense vibronic transitions (grey shadow in Figure 48a and Figure 48b) correspond to the same vibrational modes, i.e. 3 and 9. However, the spectra generated with PCM model are sharper due to the lower intensity of these medium vibronic transitions compared to the most dominant (0-0' and 0'-0) ones. Then if we compare the microsolvated systems in GP (µSolv-GP) and with PCM model (µSolv-PCM),

we found that the 0-0' and 0'-0 transitions are the most dominant ones, being the remaining vibronic transitions are quite low intense, especially the ones recorded with PCM model (Figure 48a and Figure 48b).

P a g e | 117

We can conclude that considering a microsolvation of water molecules and therefore taking into account the H-bond interactions does not modify the sharper shape of the vibronic absorption and emission spectra. However by comparing the µSolv-GP and µSolv-PCM spectra, we clearly observe that the presence of a large number of medium-intense vibronic transitions in GP broads the vibronic spectra (Figure 48a and Figure 48b).

Regarding the medium-intense band (beige shadow in Figure 48a and Figure 48b), we observe that the band obtained in not microsolvated GP is more intense than the one in not microsolvated PCM. This band presents similar amount of vibronic transitions in GP and in PCM (black and red sticks in Figure 48a and Figure 48b). While, the relative intensity of the vibronic transitions obtained in PCM is quite low compared to the most dominant 0-0' or 0'-0 transitions leading to a less intense band. In general when the microsolvation effect is included, the number of vibronic transitions reduces, but their relative intensity increases keeping the overall intensity of the medium-band similar to the one obtained in the not microsolvated systems.

We can conclude that considering a microsolvation of water molecules and therefore taking into account the H-bond interactions does not modify the sharper shape of the vibronic absorption and emission spectra. However, a significant difference in the relative intensities of the vibronic bands is observed in the computed spectra in GP or with the PCM model.

Conclusion

In this work, we have tested different solvation models to compute vertical transition energies and to simulate the absorption and emission spectra of oxyluciferin natural forms and their analogues. In addition, we have compared the wavelength maxima of the simulated spectra to the ones measured experimentally. 16,[START_REF] Garcia-Iriepa | Simulation and Analysis of the Spectroscopic Properties of Oxyluciferin and Its Analogues in Water[END_REF] The conclusion to be drawn from this work is that the nature of the solvation model used to compute the vertical transition energy of a chemical compound could be crucial to reproduce a result in agreement with the experiment. We have demonstrated the importance of the solvation effect on the vertical transition energies of the oxyluciferin natural forms and their analogues. Various behaviours have been found depending on the charge transfer (CT)

character of the systems under study. For example, the implicit SS-PCM formalism is the most suitable to compute the vertical transition energies of chemical compounds with a large CT character while, the implicit LR-PCM formalism is appropriate for compounds characterised by a small (negligible) CT character, in line with the observations of previous studies. 146,147 However, P a g e | 118

the absolute vertical transition energies computed with the implicit solvent models still do not match the experimental results, unlike the transition energies computed with the explicit solvent model which fit the experimental results.

Finally apart from these absorption and emission energies, we have also studied the solvation effect on the vibronic spectra of the phenolate-keto oxyluciferin form and explored the spectral shapes. It has been found that the relative intensities of the bands are different in the spectra simulated in gas phase or with the implicit model. However, applying a microsolvation model (considering explicitly a small sphere of solvent molecules around the solute) to include the H-bonding interactions, does not have a significant impact on the overall spectral shape.

Context

As we have seen previously (Chapter 6), the oxyluciferin natural forms can be interconverted by inter-exchange reactions such as keto-enol tautomerization or acid-base reactions (Figure 2). Hence, the photochemical properties of the oxyluciferin forms have been found to be dependent on pH value. 5 The change of the pH value can also modify the oxyluciferin surrounding, for instance the protonation state of both the protein active site and the adenosine monophosphate (AMP). The change of the protonation state of the AMP inside the protein cavity has been already investigated by using the fragment molecular orbital calculations. 27 According to the results of this investigation, the change of AMP protonation state can be responsible for the phenolate-keto polarisation by the protein active site. 27 The interplay between the protonation states of oxyluciferin-AMP complex inside the protein cavity has been studied by C. García-Iriepa and I. Navizet 28 using QM/MM methods. They have been found that the AMP deprotonation leads to a blue-shift of the emission spectrum. [START_REF] Garcia-Iriepa | Effect of Protein Conformation and AMP Protonation State on Fireflies' Bioluminescent Emission[END_REF] The AMP (chemical formula: C10H14N5O7P) is located inside the protein cavity and it is close to the thiazolone moiety of oxyluciferin structure. AMP can exist in two protonation states, mono-deprotonated (AMPH -

) and doubly-deprotonated (AMP 2-), and hence it can form different H-bonding with the neutral oxyluciferin analogue and its deprotonated form. Then, the change of its protonation state could affect the chemical form of the oxyluciferin and therefore modify the absorption energies.

The purpose of this study is to demonstrate whether the presence of the AMP at different pH values modifies the absorption spectra of the oxyluciferin analogues in water. For this aim, we evaluated the effect of the AMP protonation state on the absorption spectra of three oxyluciferin analogues (phenol-cycle, phenol-OMe and OMe-enol forms) in water at P a g e | 131 different pH values. Each oxyluciferin analogue form can be deprotonated only in a single position (Figure 49). Indeed, AMP is a good candidate to accept or donate protons. To achieve this work, our collaborators have measured experimentally and simulated theoretically the absorption spectra of the oxyluciferin analogues-AMP complexes (Figure 49) in water at different pH values (pH 5 to 11). The experimental spectra have been measured by steady-state spectroscopy and the theoretical ones have been simulated using a new protocol (recently proposed by Pr. N. Ferré group) combining Constant-pH molecular dynamics and hybrid QM/MM calculations [START_REF] Chen | Systematic theoretical investigation on the light emitter of firefly[END_REF] (CpHMD-then-QM/MM). 148,149 For our part, we have computed the vibronic absorption spectra of the oxyluciferin analogues-AMP complexes to include the vibronic contributions of the chromophores to the simulated CpHMD-then-QM/MM absorption spectra. We have applied the FCHT approximation to investigate the vibronic contributions in gas phase (GP) and in water using an implicit PCM model.

The results of this study show that the presence of AMP does not significantly shift the absorption spectra of the oxyluciferin analogues in water, whatever the pH value. [START_REF] Chen | Systematic theoretical investigation on the light emitter of firefly[END_REF] CpHMD is a method that permits to sample the conformation space while simultaneously sampling the protonation states space.

Results and discussion

Vibronic absorption spectra of oxyluciferin analogues in gas phase (GP) and in water

In order to evaluate the intramolecular vibronic effects of the oxyluciferin analogue forms (Figure 49), we have computed their vibronic absorption spectra in gas phase (GP) and in water (PCM) [START_REF] Naumov | Spectral-Structural Effects of the Keto-Enol-Enolate and Phenol Phenolate Equilibria of Oxyluciferin[END_REF] using the FCHT approximation within the Adiabatic Hessian (AH) model (Figure 50 and Figure 52). The experimental and CpHMD-then-QMMM results obtained in water and the wavelengths of the absorption maxima of the absorption spectra computed considering or not the vibrational information in GP are regrouped in Figure 51.

The obtained results show that the wavelength of the vibronic spectrum maximum (λmax,vib) of each analogue form can be associated with the corresponding 0-0' transition energy (Figure 50 and Figure 51). The absorption vibronic wavelength of the deprotonated OMe-enol form is strongly redshifted with respect to the protonated form (by 268 nm red-shift), much larger than the one obtained in water, either experimentally measured or theoretically calculated (Figure 50 and Figure 51). The red-shift observed in the case of the deprotonated phenol-OMe is less significant than the one of protonated form (135 nm), still much larger than the one observed in water. Furthermore, the absorption wavelength of the deprotonated phenol-OMe analogue form is red-shifted by 65 nm with respect to the protonated form. The protonated and deprotonated phenol-OMe forms are red-shifted by 16 nm and 144 nm with respect to the ones obtained with CpHMD-then-QMMM in water, respectively (Figure 51). For the protonated and deprotonated OMe-enol analogues, their absorption wavelengths are more red-shifted: 24 nm and 268 nm with respect to those obtained with CpHMD-then-QMMM in water, respectively (Figure 51). [START_REF] Naumov | Spectral-Structural Effects of the Keto-Enol-Enolate and Phenol Phenolate Equilibria of Oxyluciferin[END_REF] It should be noted that all the results of the vibronic spectra simulation obtained using PCM model are not part of the published work. The homogeneous broadening effect can be estimated as the energy difference between the TDDFT vertical excitation (λmax,abs) and the 0-0' band maximum (0-0' λvib) computed in gas phase: about 0.31 eV (36 nm) for the protonated phenol-OMe form and 0.17 eV (38 nm) for its deprotonated form (Figure 51). Assuming that the differential effect of the vibronic coupling between these protonated and deprotonated forms (about 2 nm or 0.14 eV) can be introduced to the CpHMD-then-QM/MM results in water (Figure 51), the total protonated to deprotonated shift of calculated λmax,abs with CpHMD-then-QM/MM increases to 9 nm (0.2 eV) instead of 7 nm (0.06 eV), still in disagreement with the experimental value of 40 nm (0.33 eV) (Figure 51). If we apply the same reasoning to the OMe-enol analogue, the total protonated to deprotonated shift of computed λmax,abs with CpHMD-then-QM/MM increases to 57 nm (0.28 eV) instead of 24 nm (0.19 eV), being in line with the experimental value of 44 nm (0.35 eV) (Figure 51). Furthermore, we have simulated the vibronic spectra of oxyluciferin analogues (phenol-OMe and OMe-enol as well as their deprotonated forms) in PCM (Figure 52). We observed that the predicted 0-0' vibronic transitions of the deprotonated phenol-OMe and the protonated OMe-enol forms are the most intense one. For the protonated phenol-OMe and the deprotonated OMe-enol forms, the predicted 0-0' vibronic transitions are no longer the dominant transitions, but still close to the wavelength spectrum maximum (Figure 52). The wavelengths of the vibronic spectrum maximum of the deprotonated form of phenol-OMe and OMe-enol analogues are redshifted, by 112 and 131 nm respectively, with respect to the protonated form (Figure 52). These redshifts are less significant with respect to the ones observed in gas phase, but still larger than the one obtained experimentally or with CpHMDthen-QMMM method (Figure 51 and Figure 52). Furthermore, the protonated and deprotonated phenol-OMe forms are red-shifted by 68 nm and 173 nm with respect to the ones obtained with CpHMD-then-QMMM in water, respectively (Figure 51 and Figure 52). For the protonated and deprotonated OMe-enol analogues, their absorption wavelengths are more red-shifted: 79 nm and 186 nm with respect to those obtained with CpHMD-then-QMMM in water, respectively (Figure 51 and Figure 52). It should be noted that before generating the vibronic spectra (Figure 52) in implicit solvent (PCM) within TDDFT, the geometry optimisations of the ground (GS) and the excited (ES) states as long as the frequencies calculations were performed with solvent in equilibrium (fully relaxed solvent fields), as Gaussian software does not give the possibility to optimise the ES with nonequilibrium solvent. Accordingly, in the ES optimisation for the calculations of the 0-0' energies and the vibrational overlaps, the solvent polarisation is determined with respect to the solute electronic density at the ES and hence it can be very different from that in the GS. Thus, the vibrational overlaps of the solvent molecules in the GS and the ES need to be included in the 0-0' vibrational overlap which is not practically feasible. To overcome this problematic, one can separate the evaluations of the solvent and the intramolecular vibronic effects of the solute as follows:

𝐸 0-0' 𝑒𝑥 = ∆𝐸 water 𝑣𝑒𝑟𝑡𝑖𝑐𝑎𝑙 -∆𝐸 GP 𝑟𝑒𝑜𝑟𝑔,𝑒𝑥 + ∆𝐸 GP 𝑧𝑝𝑒 Eq. 89

With ∆𝐸 GP 𝑟𝑒𝑜𝑟𝑔.𝑒𝑥 = 𝐸 GP 𝑒𝑥 (𝑅 𝑣𝑒𝑟𝑡𝑖𝑐𝑎𝑙 ) -𝐸 GP 𝑒𝑥 (𝑅 𝑚𝑖𝑛𝑖𝑚𝑢𝑚 ) Eq. 90

The first term (∆𝐸 water 𝑣𝑒𝑟𝑡𝑖𝑐𝑎𝑙 ) represents the vertical excitation energy in water. The second term (∆𝐸 GP 𝑟𝑒𝑜𝑟𝑔,𝑒𝑥 ) represents the intramolecular vibrational reorganisation energy in gas phase (GP).

This term gives the energy difference between the ES energies (𝐸 GP 𝑒𝑥 ) at the vertical geometry (𝑅 𝑣𝑒𝑟𝑡𝑖𝑐𝑎𝑙 ) and the minimum one (𝑅 𝑚𝑖𝑛𝑖𝑚𝑢𝑚 ) in GP. The last term (∆𝐸 GP 𝑧𝑝𝑒 ) is the zero point P a g e | 136

energies corrections of the GS and the ES computed in GP. The first term represents the FC principle of the solvents in the classical approximation, while the last two terms represent the contribution of the intramolecular structural changes of the solute.

Applying the formula of Eq. 89 for the calculation of the vibrational 0-0' excitation energy for phenol-OMe analogue, we obtain 3.02 eV (411 nm) and 2.29 eV (540 nm) for its protonated and deprotonated forms, respectively. The homogeneous broadening effect can be estimated as the energy difference between the TDDFT vertical excitation energy (λmax,abs) computed in water and the 0-0' vibrational energy (𝐸 0-0' 𝑒𝑥 ): about 0.06 eV (9 nm) for the protonated phenol-OMe form and 0.07 eV (15 nm) for its deprotonated form (Figure 51 and Figure 52). Assuming that the differential effect of the vibronic coupling between these protonated and deprotonated forms (about 6 nm or 0.01 eV) can be introduced to the CpHMDthen-QM/MM results in water (Figure 51), the total protonated to deprotonated shift of calculated λmax,abs with CpHMD-then-QM/MM increases to 13 nm (0.07 eV) instead of 7 nm (0.06 eV), still also far from the experimental value of 40 nm (0.33 eV) (Figure 51). In the case of OMeenol analogue, the 0-0' vibrational energies (𝐸 0-0' 𝑒𝑥 ) of the protonated and the deprotonated forms are 2.94 eV (422 nm) and 2.25 eV (551 nm), respectively. If we apply the same reasoning to the OMe-enol analogue, the total protonated to deprotonated shift of computed λmax,abs with CpHMD-then-QM/MM increases to 27 nm (0.17 eV) instead of 24 nm (0.19 eV), still far from the experimental value of 44 nm (0.35 eV) (Figure 51).

In summary, the absorption spectral shapes of the oxyluciferin analogues simulated with CpHMD-then-QM/MM approach in water, at different pH values are in line with the experimental ones. 150 However, small shifts of the wavelengths of the spectra maximum have been observed (Figure 51). In addition, the experimental protonated-deprotonated energy gap of the phenol-cycle (93 nm) is close to the one computed using the CpHMD-then-QM/MM approach (85 nm) (Figure 51). On the other hand, the protonated-deprotonated energy gaps of phenol-OMe and OMe-enol forms are largely underestimated (7 nm and 24 nm) using the CpHMD-then-QM/MM approach with respect to the experiment ones (40 nm and 44 nm), respectively. Including the vibronic contributions of the analogue structures to the simulated CpHMD-then-QM/MM spectra slightly improve the protonated-deprotonated energy gaps, but still in disagreement with the experimental one. Applying this correction to the CpHMD-then-QMMM spectra would only modify the peak position, not their spectral shapes.

P a g e | 137

AMP effect on the vibronic absorption spectra in water

We have studied the AMP (mono-and doubly-deprotonated) effect on the simulated absorption spectra of oxyluciferin analogues (phenol-OMe and OMe-enol as well as their deprotonated forms) by including or not the vibrational information in implicit water (PCM) [START_REF] Branchini | Yellow-green and red firefly bioluminescence from 5,5-dimethyloxyluciferin[END_REF] .

In order to simulate the vibronic spectra we applied the Vertical Gradient (VG) model [START_REF] Rebarz | Deciphering the protonation and tautomeric equilibria of firefly oxyluciferin by molecular engineering and multivariate curve resolution[END_REF] instead the Adiabatic (AH) one as the relative positions of AMP and oxyluciferin analogue are not the same in the ground and excited states (AH calculations are infeasible due to the flexible complexes including the analogue form and AMP).

Before analysing the AMP effect on the simulated absorption spectra, we compared the computed spectra of phenol-OMe and OMe-enol analogues in PCM obtained using AH and VG models (Figure 52, Figure 53a and In general, the spectra simulated considering AMP (both mono-and doubly-deprotonated) show some vibrational transitions close to the 0-0' one. However, when the spectra were simulated with the analogues alone, these transitions have not been observed. Moreover, if we compare the spectra of the [START_REF] Branchini | Yellow-green and red firefly bioluminescence from 5,5-dimethyloxyluciferin[END_REF] It should be noted that all the results of the vibronic spectra simulation obtained using PCM and Vertical Gradient (VG) model are not part of the published work. [START_REF] Rebarz | Deciphering the protonation and tautomeric equilibria of firefly oxyluciferin by molecular engineering and multivariate curve resolution[END_REF] In the VG model, the potential energy surfaces (PES) of both GS and ES are assumed to be similar, so that only the information at GS were used and the normal mode and frequencies at the ES are ignored. Furthermore, the relative intensities of most intense vibronic transitions 0-0' band of the complexes analogue-AMP are quite similar to the one of the analogues alone, except for the complex OMe-enol-doubly-deprotonated AMP for which it is half of OMe-enol (Table 4).

Moreover, we found the same trend even in the less intense vibronic transitions (see the zoomed parts). The presence of AMP in its two protonation states (mono-or doublydeprotonated AMP) red-shifts the vibronic absorption spectra. A red-shift of 100 nm have been shown in the presence of doubly-deprotonated AMP in both neutral analogues (Figure 53, Figure

P a g e | 139

A-2 and Table 4). In contrast, a less important red-shift (17 nm for phenol-OMe and 48 nm OMeenol) have been obtained in both natural analogues considering mono-deprotonated AMP. In the cases of the deprotonated forms analogues considering mono-or doubly-deprotonated AMP, a smaller red-shift (12-33 nm) with respect to the deprotonated analogues alone have been observed (Figure 53, Figure A-2 and Table 4). Furthermore, the spectra simulated considering doubly-deprotonated AMP are more red-shifted than those obtained with monodeprotonated AMP with respect to the analogues alone, except in case of the OMe-enol deprotonated form in which with mono-deprotonated AMP are more red-shifted than with doubly-deprotonated AMP (Figure 53, Figure A-2 and Table 4).

Table 4 Wavelengths values corresponding to both maxima of the absorption (λmax,abs) and the vibronic (0-0' λVib) spectra simulated in water (PCM) using VG model. Vibronic transition 0-0' band intensities are given in molar absorption coefficient (in dm 3 .mol -1 .cm -1 ).

Moreover, the absorption wavelengths (λmax,abs) simulated considering AMP (both mono-and doubly-deprotonated AMP) are shifted (Table 4). Starting with the protonated analogues forms, a red-shift of 92-105 nm is observed in presence of doubly-deprotonated AMP. However, for the complex phenol-OMe-mono-deprotonated AMP the red-shift is less important (30 nm) and negligible in the complex OMe-enol-mono-deprotonated AMP (3 nm).

For the deprotonated analogues forms, the red-shift of the wavelengths is quite small (2-5 nm)

when considering the doubly-deprotonated AMP. In contrast, the absorption wavelength maxima are blue-shifted (9-41 nm) in presence of mono-deprotonated AMP (Table 4).

Absorption spectra of oxyluciferin analogue-AMP complex in water

The experimental and CpHMD-then-QMMM absorption spectra of the oxyluciferin analogue-AMP complexes in water at different pH values obtained by our collaborators are presented in Figure 54. The complexation of AMP with the oxyluciferin analogue form in water is supposed to modify the absorption spectra, at least because of its protonation state change P a g e | 140

(-1 or -2 electric charge of AMP) in the considered pH values (pH 5 to 11). However, comparing the results presented in Figure 51 and Figure 54, we observe that the wavelength maxima of the experimental spectra of the complexes analogue-AMP are similar to the ones obtained for the analogue alone. The wavelength maxima of the CpHMD-then-QMMM computed spectra do not significantly change when AMP is included. The simulated absorption spectra of the complex phenol-cycle-AMP in the pH range are in good agreement with the experiment. However, a slight blue-shift of the QM/MM isobestic point is observed with respect to the experimental one (Figure 54a). The computed wavelengths of the absorption spectra maxima of the phenol-OMe analogue show a limited effect of AMP, through a 16 nm red-shift for the protonated (pH 5) form and 3 blue-shift for the deprotonated (pH 9) with respect to the experiment Figure 54b).

Moreover, the presence of AMP induces a slight shift of the computed wavelengths of the absorption maxima of the OMe-enol analogue with respect to the experiment (Figure 54c). The protonated-deprotonated energy gap computed for the complex analogue-AMP are similar to the one computed for the analogue alone. The protonated-deprotonated energy gaps of the phenol-OMe-AMP and OMe-enol-AMP complexes using the CpHMD-then-QMMM method are only half the experimental one (Figure 54).

P a g e | 141

Apart the slight shift of the spectra mentioned above, the complexation of AMP with the oxyluciferin analogue form in water do not significantly modify the absorption spectra, either experimentally or theoretically. So, we can state that the computed and measured absorption spectra in water are independent of AMP.

Conclusion

In this work, the absorption spectra pH-dependency for three oxyluciferin analogue forms (phenol-cycle, phenol-OMe and OMe-enol) in water and in presence or not of AMP have been investigated using experimental and theoretical protocols. It has been shown that the absorption spectra simulated in water at different pH values using a new CpHMD-then-QM/MM approach are generally in agreement with the experiment. However, the protonateddeprotonated energy gap of the analogue forms are underestimated using CpHMD-then-QM/MM approach. Furthermore, despite the fact that the vibronic contributions of the analogue structure are introduced to the simulated CpHMD-then-QM/MM spectra, the protonated-deprotonated energy gap is still in disagreement with the experimental one.

Regarding the AMP effect on the simulated absorption spectra including or not the vibronic information, it has been found that AMP does not change the spectral shapes of the computed spectra in implicit water (PCM) but it induces a shift of their picks. However, the investigations on the effect of the presence of the AMP conducted by our collaborators on the absorption spectra of three oxyluciferin analogues (phenol-cycle, phenol-OMe and OMe-enol) in water at different pH values have been shown that the complexation of AMP with the oxyluciferin analogue form does not modify significantly the absorption spectra, either experimentally or theoretically, unlike the PCM results. This disagreement can arise from the explicit interactions between the solute and the water molecules considered when using CpHMD-then-QMMM method, which cannot be considered in isotropic solvent PCM models. Firefly's luciferase and luciferin are widely used for bioluminescence imaging, [START_REF] Woodroofe | Novel Heterocyclic Analogues of Firefly Luciferin[END_REF][START_REF] Greer | Imaging of light emission from the expression of luciferases P a g e | 203 in living cells and organisms: a review[END_REF][START_REF] Kuchimaru | A luciferin analogue generating nearinfrared bioluminescence achieves highly sensitive deep-tissue imaging[END_REF] for instance for detecting tumours in vivo. The efficiency of this technique depends on the brightness of the signal, that is, the emitted bioluminescent light. However, the emitted light from luciferin-luciferase system is yellow-green (around 550-560 nm) in colour which can be readily absorbed by the haemoglobin and hence is not able penetrate biological tissues. To overcome this obstacle and improve the signal detection, several experiments have been done to design a system emitting light colour ranged between red and near-IR spectral window (> 600 nm) as it has been shown that wavelengths of this spectral range are strong enough to penetrate deeper living tissues. [32][33][START_REF] Weissleder | A clearer vision for in vivo imaging[END_REF] Among the possibilities to tune the light colour of firefly bioluminescence from yellow to near-IR spectral window we can mention the chemical modification of the light emitter structure [START_REF] Woodroofe | Novel Heterocyclic Analogues of Firefly Luciferin[END_REF][START_REF] Ioka | Synthesis of Firefly Luciferin Analogues and Evaluation of the Luminescent Properties[END_REF][START_REF] Ikeda | An allylated firefly luciferin analogue with luciferase specific response in living cells[END_REF] and the mutation of the luciferase [START_REF] Mofford | Aminoluciferins Extend Firefly Luciferase Bioluminescence into the Near-Infrared and Can Be Preferred Substrates over D-Luciferin[END_REF][START_REF] Wang | Impact of site-directed mutant luciferase on quantitative Green and Orange/Red Emission Intensities in Firefly Bioluminescence[END_REF][START_REF] Jones | Orthogonal Luciferase-Luciferin Pairs for Bioluminescence Imaging[END_REF][START_REF] Branchini | Red-emitting luciferases for bioluminescence reporter and imaging applications[END_REF][START_REF] Branchini | Chemically Modified Firefly Luciferase Is an Efficient Source of Near-Infrared Light[END_REF] or both of them together. Despite the efforts made to clarify the colour tuning phenomenon, many details are still unsolved. [START_REF] Navizet | Color-Tuning Mechanism of Firefly Investigated by Multi Configurational Perturbation Method[END_REF][START_REF] Hosseinkhani | Molecular enigma of multicolor bioluminescence of firefly luciferase[END_REF][START_REF] Iwano | Development of simple firefly luciferin analogs emitting blue, green, red, and near-infrared biological window light[END_REF] In the case of the chemical modification of the light emitter structure, various synthetic analogues have been designed to obtain a red/near-IR emission colour. However, the chemical structure modifications are not an easy task and, in some cases, the emission wavelength of analogues can be observed in blue and in red colours or even outside of the colour spectral window, depending on the substitution made in the light emitter structure moieties, benzothiazole [START_REF] Woodroofe | Novel Heterocyclic Analogues of Firefly Luciferin[END_REF][START_REF] Ikeda | An allylated firefly luciferin analogue with luciferase specific response in living cells[END_REF][START_REF] Zhang | Pyridone Luciferins and Mutant Luciferases for Bioluminescence Imaging[END_REF] and thiazolone [START_REF] Ioka | Synthesis of Firefly Luciferin Analogues and Evaluation of the Luminescent Properties[END_REF][START_REF] Mccutcheon | Expedient synthesis of electronically modified luciferins for bioluminescence imaging[END_REF][START_REF] Jathoul | A dual-color far-red to nearinfrared firefly luciferin analogue designed for multiparametric bioluminescence imaging[END_REF] or in the π-backbone. [START_REF] Iwano | Development of simple firefly luciferin analogs emitting blue, green, red, and near-infrared biological window light[END_REF][START_REF] Jathoul | A dual-color far-red to nearinfrared firefly luciferin analogue designed for multiparametric bioluminescence imaging[END_REF] The purpose of this work is to investigate the influence of structural modifications within the firefly light emitter structure, oxyluciferin, on the colour tuning of bioluminescence. For this P a g e | 155

aim, we studied the natural firefly oxyluciferin (in its phenolate-keto [START_REF] Naumov | Spectral-Structural Effects of the Keto-Enol-Enolate and Phenol Phenolate Equilibria of Oxyluciferin[END_REF] chemical form, hereafter donated Natural_Oxy) and three of its synthetic analogues named benzothiophene_Oxy, dihydropyrrolone_Oxy and allylbenzothiazole_Oxy (Figure 55). In the two first analogues, the structural modifications are located in the benzothiazole moiety (nitrogen atom replaced by a methine) and in the thiazolone moiety (sulfur atom replaced by a methylene), respectively, exhibiting blue-shifts of the emission wavelengths with respect to the natural bioluminescence.

The last analogue is synthesised by introducing an allyl group in the benzothiazole moiety,

showing a red-shift of the emission wavelength (Figure 55). To realise this work, we performed a classical MD simulation for each compound considering explicitly the water molecules and the protein environment in order to obtain a statistical number of snapshots. Afterward, we carried out the QM/MM calculations to compute the vertical transition energies and then simulate the absorption and the emission spectra (see computational details in Chapter 5). Moreover, we computed the vertical transition energies not only considering the protein environment but also in gas phase (GP) and in implicit water (PCM) at the QM level of theory for comparison. In addition, we analysed the substituents and solvent effects on the chemical physical properties: (i) the vertical transition energies (Te) and oscillator strength (f), (ii) the charge transfer (CT) character and the natural transition orbitals (NTOs), and

(iii) the energies of the HOMO-LUMO orbitals involved in the vertical transitions. Then, we compared the predicted results to the ones measured experimentally. [START_REF] Woodroofe | Novel Heterocyclic Analogues of Firefly Luciferin[END_REF][START_REF] Ioka | Synthesis of Firefly Luciferin Analogues and Evaluation of the Luminescent Properties[END_REF][START_REF] Ikeda | An allylated firefly luciferin analogue with luciferase specific response in living cells[END_REF] In this work, we demonstrated that the experimental bioluminescence spectra can be reproduced by using the QM/MM methods including explicitly the protein environment.

Furthermore, the predicted results show that the structural modification within the oxyluciferin structure yields changes of its electronic distribution and HOMO-LUMO energies, leading to [START_REF] Naumov | Spectral-Structural Effects of the Keto-Enol-Enolate and Phenol Phenolate Equilibria of Oxyluciferin[END_REF] The choice of using the phenolate-keto form of oxyluciferin is based on assumption that it is the most probable chemical form responsible of the emission in fireflies bioluminescence. [START_REF] Chen | Systematic theoretical investigation on the light emitter of firefly[END_REF][START_REF] Naumov | Spectral-Structural Effects of the Keto-Enol-Enolate and Phenol Phenolate Equilibria of Oxyluciferin[END_REF][START_REF] Branchini | Yellow-green and red firefly bioluminescence from 5,5-dimethyloxyluciferin[END_REF][START_REF] Rebarz | Deciphering the protonation and tautomeric equilibria of firefly oxyluciferin by molecular engineering and multivariate curve resolution[END_REF][START_REF] Naumov | Spectral-Structural Effects of the Keto-Enol-Enolate and Phenol Phenolate Equilibria of Oxyluciferin[END_REF] P a g e | 156

different light emission colours. Finally, thanks to this study, we proposed a novel oxyluciferin analogue which could be a good candidate to emit light in the near-IR spectral window.

Results and discussion

Compounds in gas phase (GP) and implicit water (PCM)

The computed vertical energies of the S0→S1 (absorption) and S1→S0 (emission) transitions of the Natural_Oxy and its three analogues in gas phase (GP) and in implicit water (PCM) are regrouped in Table 5. By analysing these computed results in GP and in PCM, we found that the absorption and emission transition energies computed in PCM were shifted to lower (around 0.07-0.11 eV) and higher (around 0.04-0.10 eV) energies, respectively, with respect to the ones obtained in GP. Generally, the differences between GP and PCM transition energy values can be due to the solute-solvent interactions which can sharply stabilise/destabilise the ground (GS) or excited (ES) states. In addition, the oscillator strengths (f) values of the absorption and emission transitions increase substantially when considering the implicit water. To rationalise this finding, we calculated the square of the transition dipole moments (µge) 2 between the GS and ES for the absorption and the emission transitions, respectively, in the GP and in PCM (Table A-2). We observed that the (µge) 2 values increase with the increasing of the oscillator strengths (f) values, as the f magnitude for an electronic transition is proportional to the square of µge.

We also computed transition energies of all the compounds in PCM using DiMethylSulfoxyde (DMSO) as solvent and compared them to those computed in water solvent in order to investigate solvent polarity effect on the predicted results, as the dielectric constants (ε) of the water and DMSO solvents within PCM model are 78.355 and 46.826, respectively.

Comparing the predicted results in DMSO (PCM) and in water (PCM), we found that the computed transition energies, oscillator strengths and the maximum absorption and emission wavelengths in DMSO are quite similar to the ones obtained in water (Table 5 and Table A-3).

Therefore, the polarity of the solvents (water and DMSO) does not significantly affect the absorption and emission wavelengths, neither the f values.

Regarding the effect of the oxyluciferin structural modifications on the absorption and emission transitions, we observed that the transition energies of the Natural_Oxy and its analogues follow the same trend as the experimental bioluminescence, that is, the maximum wavelengths of the benzothiophene_Oxy and dihydropyrrolone_Oxy spectra are blue-shifted P a g e | 157

with respect to the ones of Natural_Oxy, while those of allylbenzothiazole_Oxy are red-shifted (Table 5). Unfortunately, a direct comparison between the experimental and computational results cannot be made, as the experimental spectra of the analogues under study in water are

not measured yet.

Furthermore, we computed the frontier molecular orbital energies (EHOMO and ELUMO) for the absorption and the emission transitions in GP and in PCM (Table A-4). It was found that both in absorption and emission transitions, the HOMO-LUMO energy differences (∆EHOMO-LUMO) were affected by the chemical modifications within the oxyluciferin structure, which is in agreement with the differences observed in the absorption and emission energies (Table 5). For instance, the ∆EHOMO-LUMO of the allylbenzothiazole_Oxy analogue is smaller than the one of the Natural_Oxy, which explaines the red-shift of the absorption and emission spectra (Table 5).

While, the ∆EHOMO-LUMO of the benzothiophene_Oxy and dihydropyrrolone_Oxy analogues are larger than the one of the Natural_Oxy, which is in line with the observed blue-shift of the absorption and emission wavelengths (Table 5). The values in brackets correspond to the wavelengths of the experimental bioluminescence spectra maximum. [START_REF] Woodroofe | Novel Heterocyclic Analogues of Firefly Luciferin[END_REF][START_REF] Ioka | Synthesis of Firefly Luciferin Analogues and Evaluation of the Luminescent Properties[END_REF][START_REF] Ikeda | An allylated firefly luciferin analogue with luciferase specific response in living cells[END_REF] (Table reproduced from ref. 151 ).

P a g e | 158

Moreover, we also analysed the charge transfer (CT) character of the absorption (S0→S1) and emission (S1→S0) transitions of the Natural_Oxy and its analogues (Table A -5). We observed that the CT character of the Natural_Oxy structure and its analogues are similar, being slightly smaller for dihydropyrrolone_Oxy and benzothiophene_Oxy structures (Table A -5). This indicates that the oxyluciferin structural modifications barely influence the electronic nature of the vertical transitions, that is, the charge transfer character.

Finally, we computed the natural transition orbitals (NTOs) involved in the absorption and emission transitions in GP and in PCM for the Natural_Oxy and its analogues (Figure 56 and group is slightly larger than the one located in N atom of the Natural_Oxy (dashed blue circle in Figure 56). Regarding the dihydropyrrolone_Oxy (where CH2 group replaces the S atom of the Natural_Oxy thiazolone moiety), the unoccupied and occupied NTOs electron densities located in the CH2 group are smaller than the one located in S atom of the Natural_Oxy (dashed blue circle in Figure 56) due to the larger electronegativity of the S atom. In the case of allylbenzothiazole_Oxy, although the introduced allyl group is not part of the π-conjugated system, its occupied NTO has some electronic density (dashed blue circle in Figure 56). Hence, we can conclude that the oxyluciferin structural modification slightly affects the NTO electron densities of the absorption and emission transitions, which is in line with the small variations of the charge transfer character observed above. The difference in the charge transfer characters and electronic densities of the compound structures is not so large but consistent with the trend of the computed absorption and emission wavelengths.

Compounds in protein

To investigate in more detail the effects of the oxyluciferin structural modification on the experimental colour modulation of the firefly's bioluminescence, we computed the vertical energies of the S0→S1 (absorption) and S1→S0 (emission) transitions of the Natural_Oxy and its analogues at the QM/MM level considering the protein environment (Table 5) and we simulated their absorption and emission spectra. By comparing the computed results in GP and in protein, we observed that the absorption and emission transition energies computed in protein were shifted to lower energies (around 0.01-0.08 eV for the absorption and 0.05-0.2 eV for the emission) with respect to those obtained in GP (Table 5). In addition, the oscillator strengths (f) of the absorption and emission transitions computed in protein are in the same order of magnitude as those in GP, being slightly smaller when considering the protein environment than in GP. The relative order of the simulated maximum absorption and emission wavelengths in GP and in PCM are the same as the ones computed in the protein (Table 5). Moreover, the natural transition orbitals (NTOs) of the Natural_Oxy and its synthetic analogues inside the protein are quite similar to those computed in GP and in PCM (Figure A-3).

For the purpose of determining the structural modification effect on the firefly's bioluminescence, we compared the emission spectra of the Natural_Oxy and its analogues simulated in protein to the ones measured experimentally. [START_REF] Woodroofe | Novel Heterocyclic Analogues of Firefly Luciferin[END_REF][START_REF] Ioka | Synthesis of Firefly Luciferin Analogues and Evaluation of the Luminescent Properties[END_REF][START_REF] Ikeda | An allylated firefly luciferin analogue with luciferase specific response in living cells[END_REF] The spectral shapes of the simulated spectra are similar to the experimental ones (Figure 57). In addition, the simulated maximum emission wavelengths are only shifted by less than 10 nm with respect to the experiment. Moreover, the maximum emission wavelengths of the simulated spectra show qualitatively the same trend as the experimental ones (Figure 57). For instance, for P a g e | 160

benzothiophene_Oxy and dihydropyrrolone_Oxy analogues, the maximum emission wavelengths are blue-shifted theoretically by 40 nm and 11 nm and experimentally by 37 nm and 13 nm, respectively, with respect to the Natural_Oxy. For the allylbenzothiazole_Oxy analogue, the simulated and experimental emission spectra are red-shifted by 46 nm and 45 nm, respectively, with respect to the Natural_Oxy. However, the simulated allylbenzothiazole_Oxy spectrum is broader than the other spectra, either theoretically or experimentally obtained (Figure 57). The origin of this spectral broadening can be due to the large fluctuation (between It should be noted that the relative order of the simulated absorption spectra of the Natural_Oxy and its analogues in protein (Figure A-5) is the same as the ones in the simulated and experimental emission spectra. Unfortunately, a direct comparison between the experimental and computational results cannot be made, as the experimental absorption spectra of the analogues under study are not available.

P a g e | 161

Furthermore, we investigated H-bonding between the Natural_Oxy (and its analogues) and the protein active site environment to better understand the influence of the oxyluciferin structural modifications. To determine the effect of the oxyluciferin structural modification on the H-bonding networks, we calculated a statistical average number of H-bonding between Natural_Oxy (and its analogues), water molecules, AMPH and the protein active site along the MD simulations of the system in the excited state (ES) (Figure 58a). A representative MD snapshot of each compound under study inside the protein active site showing the H-bonding is presented in Figure 58c. If we compare the H-bonding formed in the case of Natural_Oxy to that of benzothiophene_Oxy analogue, we observe that the Natural_Oxy H-bonding formed between the N2 atom (see Figure 58b for atom numbering) of the benzothiazole moiety and the SER347 and GLY246 amino acids does not exist anymore due to the structural modification (Figure 58c).

Moreover, we found that benzothiophene_Oxy analogue gets closer to AMPH structure leading to more H-bonding than in the case of the Natural_Oxy (Figure 58c and Figure A-6a). In addition, the Natural_Oxy H-bonding formed between the oxygen atom (O10) of the phenol group and ARG337 no longer exists in the case of benzothiophene_Oxy analogue (Figure 58a andc), which can be due to the significant flexibility of the ARG amino acid side chain as already observed in previous work with Natural_Oxy. [START_REF] Garcia-Iriepa | Effect of Protein Conformation and AMP Protonation State on Fireflies' Bioluminescent Emission[END_REF] To confirm this hypothesis, we performed a second MD simulation for the benzothiophene_Oxy analogue (Figure 59c). This second MD simulation shows that ARG337 amino acid comes close to the oxygen atom (O10) of the phenol group leading to the H-bonding. Hence, the structural modification can be excluded as reason for this finding.

P a g e | 162 simulations (Figure 59). We found that the measured distance between dihydropyrrolone_Oxy atoms and the amino acids or AMPH atoms are located almost at the same distances as the measured one for the Natural_Oxy (Figure 59). In addition, the calculated statistical average number of H-bonding formed during the MD simulations of both dihydropyrrolone_Oxy and Natural_Oxy are quite similar (Figure 58a). Hence, the oxyluciferin structural modification does not significantly affect the interaction between the dihydropyrrolone_Oxy structure and its surrounding. This means that the origin of the emission spectra blue shifts does not lie on different formed H-bonding.

The H-bondings formed between the allylbenzothiazole_Oxy analogue and either amino acids (SER347, ARG337, ARG218 and GLY246) or AMPH during the MD simulation are similar to those formed with the Natural_Oxy. However, the H-bond network formed between the oxygen atom (O10) of the phenol group and SER314 amino acid (close to the allyl group) no longer exists in the case of the allylbenzothiazole_Oxy analogue (Figure 58c). The measured distance between O10 atom of allylbenzothiazole_Oxy and the SER314 amino acid is around 7.5 Å, while the measured one for the Natural_Oxy is around 5.5 Å (Figure 59a). In addition, the calculated statistical average number of H-bondings formed between the O10 atom of allylbenzothiazole_Oxy analogue and the surrounding water molecules is half of that for Natural_Oxy (Figure 58a The conclusion to be drawn from this part is that the oxyluciferin structural modification slightly influences on the one hand, the distribution of the electronic density, and on the other hand, the H-bonding between the compounds and its surrounding: the protein environment, water molecules and AMPH. However, the formed H-bonding networks during MD simulations did not cause the shifts of the emission spectra, as is the case of dihydropyrrolone_Oxy spectrum blue-shift. Finally, the emission spectra of the Natural_Oxy and its analogues simulated considering the protein environment by using the QM/MM approach are in quite good agreement with the experimental bioluminescence ones.

Suggestion of novel oxyluciferin analogue

Based on the work done so far, we may propose a novel oxyluciferin analogue in two conformations, Red_1_Oxy and Red_2_Oxy analogues (Table 6), which could be a promising candidate for emitting light in the red/near-IR spectral window. This analogue has been conceived by combining two modifications with the oxyluciferin structure: the first modification consists of introducing an allyl group on the phenol group of the oxyluciferin benzothiazole P a g e | 165 moiety, and the second one by extending the π-conjugated system, adding a double bond between the oxyluciferin moieties (Table 6).

With the aim of determining if this new analogue is a good candidate for emitting light in the near-IR spectral window, we computed the emission transition energies (Te) and the corresponding wavelengths (λemi) for both of its conformations, Red_1_Oxy and Red_2_Oxy, in gas phase (GP) and in implicit water (PCM), as we have previously seen that reasonable calculations in GP and in PCM can already provide a good qualitative trend of the emission spectra (Table 6). If we compare the predicted results of this novel analogue to those of the Natural_Oxy and allylbenzothiaozole_Oxy computed in GP and in PCM (Table 5 and Table 6), we observe that the maximum emission wavelengths of two conformations of the new analogue are red-shifted with respect to those of Natural_Oxy and allylbenzothiaozole_Oxy analogue. For instance, the maximum emission wavelengths of the Red_1_Oxy form are red-shifted by 54 nm (in GP) and 62 nm (in PCM) compared to the Natural_Oxy. For the Red_2_Oxy form, the maximum emission wavelengths are red-shifted by 83 nm (in GP) and 89 nm (in PCM) with respect to Natural_Oxy ones, being more red-shifted (29 nm in GP and 27 nm in PCM) than the Red_1_Oxy form (Table 5 and Table 6).

Table 6 Vertical transition energies (Te), the corresponding wavelengths (λemi) and oscillator strengths (f) of the emission transition of Red_1_Oxy and Red_2_Oxy analogues in gas phase (GP), in implicit water (PCM) and considering the protein [START_REF] Navizet | Color-Tuning Mechanism of Firefly Investigated by Multi Configurational Perturbation Method[END_REF] environment. For the compound in protein, the Te value corresponds to the maximum emission wavelength of the emission spectrum and the f values are taken from representative MD snapshots. (Table reproduced from ref. 151 ). [START_REF] Navizet | Color-Tuning Mechanism of Firefly Investigated by Multi Configurational Perturbation Method[END_REF] It should be noted that all the results obtained for Red_2_Oxy conformer in protein are not part of the published work.

P a g e | 166

Furthermore, we computed the natural transition orbitals (NTOs) involved in the emission transitions for the Red_1_Oxy and Red_2_Oxy conformers in GP and in PCM (Figure 60). By analysing the NTOs, we found that their electron density computed in GP are quite similar to the one in PCM. In addition, we observed that the electron density of the two conformers are transferred from the thiazolone moiety to the benzothiazole moiety upon emission transition.

Furthermore, the extension of the π-conjugated system of the oxyluciferin structure and introducing an allyl group on its benzothiazole moiety slightly influence the electronic density distribution which in turn affects the emission transition energies (Te), leading to the redshift of the maximum emission wavelengths (Table 5 and Table 6). Since the QM/MM calculations are needed for quantitative agreement with the experiment results, we performed MD simulation of the Red_2_Oxy conformer as it shows interesting results in GP and in PCM and computed the emission transition energies (Te) and then simulated the emission spectrum considering the protein environment at the QM/MM level. In general, considering the protein environment in the QM/MM calculations redshifts the maximum emission wavelengths. If we compare the emission energy of the Red_2_Oxy conformer in protein to the one computed in GP, we observe that the emission energy was shifted to lower energy (by 0.2 eV) when considering the protein environment (Table 6). The oscillator strength (f) of the Red_2_Oxy conformer in protein is of the same order of magnitude as the one computed in GP, being slightly smaller in the protein than in the GP. In addition, the f values of the Red_2_Oxy conformer in GP and in protein environments, are slightly larger than the ones computed for Natural_Oxy, and allylbenzothiazole_Oxy analogue (Table 5 and Table 6).

P a g e | 167

Regarding the simulated spectrum of the Red_2_Oxy conformer in protein, it presents similar spectral shape compared to those computed for the Natural_Oxy and allylbenzothiazole_Oxy (Figure 61). However, the computed maximum emission wavelength of the Red_2_Oxy conformer is red-shifted by 132 nm and 86 nm comparing to the ones of the Natural_Oxy and allylbenzothiazole_Oxy, respectively. It should be noted that the Red_2_Oxy simulated spectrum is broader than the one of Natural_Oxy (Figure 61). This may be due to the We can state that the novel suggested analogue, the Red_2_Oxy conformer, is a good promising candidate for emitting near-IR light, subject to validation by experimental measurements.

Conclusion

In this study, the vertical transition energies (Te) and the absorption and emission spectra of firefly oxyluciferin and some of its analogues that emit different light emission colours have been simulated at QM and QM/MM levels of theory. The predicted results in gas phase (GP) and in water (PCM) show that considering the solvent, the absorption and emission energies are shifted to lower (0.07-0.11 eV) and higher (0.04-0.1 eV) energies, respectively.

Moreover, the analysis of the computed absorption and emission transition energies (Te) in GP P a g e | 168

and in protein shows that considering the protein environment shifts their values to lower energies (around 0.01-0.08 eV for the absorption and 0.05-0.2 eV for the emission). The effects of the oxyluciferin structural modification on the emission values can be summed up as follows:

modifications by less electron donor groups in the thiazolone (S atom by CH2 group) moiety or in the benzothiazole (N atom by CH group) one induce a blue-shift of the absorption and emission spectra, while modification by more electron donor group on benzothiazole (H atom by CH2CHCH2 group) moiety yields a red-shift of the absorption and emission spectra. In addition, the oxyluciferin structural modification slightly influences on the one hand, the distribution of the electronic density and on the other hand, the H-bonding between the compound and its surrounding: the protein environment, water molecules and AMPH. However, the formed H-bonding networks during MD simulations did not cause the shifts of the emission spectra, as is the case of dihydropyrrolone_Oxy spectrum blue-shift. Moreover, the steric effects of induced by the allyl group affect the H-bonding formation. The emission spectra of the Natural_Oxy and its analogues simulated considering the protein environment by using the QM/MM approach have been found in quite good agreement with the experimental bioluminescence ones. Finally, the suggested novel analogue, Red_2_Oxy analogue, has been found as a promising candidate to emit near-IR light, subject to validation by experimental measurements. Thanks to the computational procedure followed during this Ph.D. Thesis it was possible to simulate reliable absorption and emission spectra of the light emitter and its analogue forms either in water solution or considering also the protein environment. This procedure consists of performing a classical molecular dynamics (MD) simulation for the system under study in order to obtain a statistical number of conformations (called snapshots in this Ph.D. Thesis) and computing the vertical transition energies of some MD snapshot using the hybrid quantum mechanics/molecular mechanics (QM/MM) methods and then simulating the absorption and emission spectra. This computational procedure gives a more realistic picture of the simulated molecular system and achieves a quite good agreement with the experimental results. The use of this procedure is a great challenge as it requires good knowledge of the system of interest and takes enormous amounts of human and computer CPU time. The most difficult task when using this computational procedure is related to the complexity to create the first model of oxyluciferin-luciferase system and to prepare a set of force field parameters for each oxyluciferin structure (or its analogues).

The fact that the structure of the firefly's light emitter exists in six different chemical forms makes the understanding of the de-excitation processes more difficult, both in theoretical or experimental researches. To deal with this complicated situation, theoretical and experimental teams work collaboratively by pooling their expertise, experience and knowledge on firefly's bioluminescent system. Through collaboration with our partners, the photochemical properties of all the possible chemical forms of the light emitter and some of their synthetic analogues have been studied in water solution to give insight into the analogue forms suitability to mimic natural light emitter forms in water. The results of this collaboration have been discussed in Chapter 6. We have demonstrated that: i) the synthetic oxyluciferin analogues P a g e | 180

represent similar H-bond networks along the MD simulation compared the ones observed for their respective natural forms, ii) the electronic nature of the vertical transitions and spectral band of the analogue forms are almost identical to that of the natural oxyluciferin forms and iii) the absorption and emission spectra of the analogue forms match the spectra of their respective natural forms. From these elements we concluded that the oxyluciferin synthetic analogues are suitable to imitate their respective natural light emitter forms in water solution.

The role of the computational solvation models on the computed vertical transition energies and hence on the simulated absorption and emission spectra of oxyluciferin natural forms and the analogues has been investigated, as the way that the solvent is treated could be crucial to reproduce a result in agreement with the experiment. The results of this study have been discussed in Chapter 7. Various behaviours have been found depending on the charge transfer (CT) character of the system under study. For instance, some solvation model formalisms are suitable for computing vertical transition energies of chemical compounds with a large CT character, while other formalisms are appropriate only for compounds characterised by a small (negligible) CT character. From this study we concluded that i) the tested implicit solvation models are not relevant to reproduce results in agreement with the experiment and ii) only the explicit model which includes the solute-solvent molecules interactions is reliable and it gives results in line with the experiment.

Another collaboration with our partners has been focused on the effect of the presence of the adenosine monophosphate (AMP) on the absorption spectra of some oxyluciferin analogues at different pH values in water solution. It has been found that the complexation of AMP with the oxyluciferin analogue form in water do not significantly modify the absorption spectra, either experimentally or theoretically. So, we concluded that the computed or the measured absorption spectra in water solution are independent of AMP. This study has been presented in Chapter 8.

The influence of structural modifications within the firefly light emitter structure, oxyluciferin, on the emission spectrum has been also studied to clarify and get insight into the firefly's bioluminescence colour tuning. A natural firefly oxyluciferin and three of its synthetic analogues that exhibit blue-and red-shifts of the emission wavelengths with respect to the natural bioluminescence have been studied in the Chapter 9. We have shown that the oxyluciferin structural modification slightly influences on the one hand, the distribution of the electronic density, and on the other hand, the H-bonding between the compound and its surrounding: the protein environment, water molecules and AMPH. Finally, we have proposed P a g e | 181

a novel oxyluciferin analogue which could be a good candidate to emit light in the near-IR spectral window.

Computational methods are essential to learn more about firefly bioluminescence system as they reveal hidden information from an experimental point of view like the electronic nature of the vertical transitions, H-bond networks between the substrate and its surrounding, etc. The choice of the computational procedure is primordial to obtain a realistic picture of the firefly's bioluminescent system and hence reproduce the experimental data. The static computational procedure (MD simulations followed by QM/MM calculations) used in this Ph.D.

Thesis to study the last step of the firefly bioluminescent reaction was seen to be effective and reliable. However, the exploration of the entire firefly bioluminescent mechanism requires the use of dynamic computational approaches. A study evaluating the performance of static and dynamic computational procedures has been started in collaboration with Roca-Sanjuán group.

This study tests distinct theoretical approaches on a small organic molecule, acrolein, in gas phase and in water (article in preparation).

Finally, series of questions still come up regarding some details in firefly's bioluminescent system for example: i) how is the light emitter positioned inside the protein active site? ii) Does it always take the same known position or can it be positioned differently?

iii) What about the AMP in these cases? vi) Does the light emitter remain or not inside the protein active site after light emission? The first three questions can be answered applying our computational procedure, it is just a matter of time and funding of course. The last question can be answered by using another computational method, Umbrella Sampling (US) method. During this Ph.D. Thesis, I have started to perform US simulation of firefly's bioluminescent system. The primary results of this simulation show that the light emitter come out of the protein active site, but in an unexpected direction (see Figure A-8). The binding energy ∆G (difference between the highest and lowest values of the potential mean of force curve) is found around 37 kcal/mol.

P a g e | 184

Table A-3 Vertical transition energies (Te), the corresponding wavelengths (λabs and λemi) and oscillator

strengths (f) of the absorption and emission of Natural_Oxy and its analogues in DMSO (PCM).

Table A-4 Frontier molecular orbital energies (EHOMO and ELUMO (a.u.)) and their energy difference (eV), computed at the absorption and emission transitions for all the structures in gas phase and in water (PCM).

Charge transfer (CT) character

To obtain CT character, we proceed as follows: we calculate, for both the ground (GS) and excited state (ES) minima, the sum of the Mulliken atomic charges included in the thiazolone (or benzothiazole) moiety of the chemical structure in the ground (S0) and in the first excited state (S1). Then, we calculate the difference between the thiazolone moiety charges in S0 and S1,

showing the through space CT character of the transition (between the thiazolone and the benzothiazole moieties) (Table A-5).

P a g e | 185

P a g e | 186

P a g e | 187

P a g e | 188

A-2 Tutorial for System preparation -Input files

In the following section, additionally to the computational details reported above (Chapter 5), we present a small tutorial explaining how the input files for the MD simulations and the QM/MM calculations have been prepared. We start by presenting input files for computing a solvation box of a given system, explaining the way in which the force field parameters are obtained, followed by performing classical molecular dynamics MD simulations.

Some tips for analysing the MD simulation will be presented. Finally, after the snapshots extraction along the MD simulation, QM/MM calculations will be carried out to compute the vertical transition energy and simulate the absorption and emission spectra.

A-2.1 Preparing the system to perform MD simulation

From the output file from Gaussian g09 "file.out" of the geometry optimisation in gas phase we extract the optimised geometry (Figure 62) and we create a coordinate file "file.xyz" keeping only the "Coordinates X Y Z" as follows: Then, with VMD we can open the "file.xyz" and save the geometry coordinates with the pdb format, given "file.pdb" as follows: When dealing with a ligand inside a protein, the classification number of its coordinates can take another classification, for example, the classification number of the oxyluciferin inside the protein cavity in our coordinate files is 544.

Once the "file.pdb" files are carefully prepared, with the Antechamber program we can create "file.prep" file for LEaP using the following command:

→ module load ambers → antechamber -i file.pdb -fi pdb -o file.prep -fo prepi P a g e | 191

The generated "file.prep" output file will look like as shown in Figure 64. The second column in this file represents the atom types (as previously fixed in file.pdb), the third column indicated the atom class (we have chosen to put numbers for further parameters and Tinker files) and the last column represents the atomic charges. In general the generated "file.prep" file, the last column is given with zero values. Thus, the user should replace them by hand with the Mulliken atomic charges of the optimised geometry that are given in the output file of Gaussian after the last SCF procedure "Population analysis using the SCF density" calculated in gas phase at the QM level. If you deal with the TDDFT calculations, that is, excited states, one needs to add the following keyword "density=current" in the Gaussian command line to get the atomic charges at the excited states. They are given at the last CI procedure "Population analysis using the CI density".

Using the Pamchk2 program which analyses "file.prep" file, we can extract the force field parameters (bond, angle …) that are needed for the molecule into force field modification file "file.frcmod".

→ parmchk2 -i file.prep -o file.frcmod -f prepi

This command reads the "file.prep" and generates the force field parameters. In some cases, it does not provide all the needed parameters. Therefore, one needs to add them manually by using the GAFF force field parameters stored in the gaff.dat file as reference (Figure 65). Once all the files (.pdb, .prep, .frcmod) are carefully prepared for the ligand, we can create the whole system, ligand with or without the protein in a solvation box, by using tLEaP program with the following command:

→ tleap leaprc

Here, "leaprc" file (Figure 66) loads the force field parameters (for the ligand and the protein), creates a pdb file with a truncated octahedron solvation box and saves the topology and coordinate files. The input file "file_MD.in" looks like:

A-2.3 Molecular dynamics simulation analysis

At the end, one needs to visualise the MD simulation with VMD to ensure that the simulation is done correctly. Besides looking at the evolution of the energies, one can compare the deformation of part of whole system during the simulation by computing some Root mean square deviation (RMSD) analysis, which provides information about the stability of the system along the MD simulation. The RMSD of the system can be computed with Cpptraj program as follows:

→ cpptraj → parm file_box.top "Reading 'oxy_box_last.top' as Amber Topology" → trajin file_MD.mdcrd "Reading 'md1_center.mdcrd' as Amber Trajectory" → rms ToFirst :1-545&!@H= first out rmsd_file.agr mass

Here, the computed mass-weighted RMSD are obtained using all non-hydrogen atoms (only the heavy atoms are considered) in the ligand (n° 544), protein residues (n° 1 to 543) and AMPH (n° 545). The structure of the first trajectory frame is taken as reference (Figure 70). In this example, the RMSD indicates that the system is stabilised after 8 ns of simulation with respect to the initial trajectory and with respect to the equilibration region. The small fluctuations in the stabilised region in the RMSD are due to the flexibility of the system during the MD simulation.

A-2.4 Hydrogen bonds (H-bonds) analysis

There are several programs to find and track the H-bonds along the MD simulation such as Cpptraj program. They permit to find not only the solute-solute H-bonds, but also solutesolvent H-bonds by using specific commands.

A-2.4.1 Solute-solute H-bonds with Cpptraj

One way of performing the solute-solute H-bonds analysis with Cpptraj program is to specify the following commands: The hbond keyword will look for H-bond donors/acceptors in region specified by: 1-545 (only the ligand "n° 544", the protein residues "from n° 1 to 543" and AMPH "n° 545" are considered, excluding the water molecules) and only look for intermolecular H-bonds within the distance cutoff = 3 Å and angle cutoff = 135°". The output file containing solute-solute H-bond averages is saved as "avg_hbond.dat" (Figure 71). In the example reported in Figure 71, we can observe from the fraction column that for 78% of the MD simulation time there is a H-bond between the oxygen atom (O10) of the ligand and the hydrogen atom (HH21) of the ARG residue. Moreover, we can also observe that 41,7% of the MD simulation time there is a H-bond between the oxygen atom (O10) of the ligand and the hydrogen atom (HE) of the ARG residue. While, we can observe that the H-bonds are short lived (0.5% of the MD simulation time) and hence, they do not yield significant stability of the The hbond keyword will search for H-bond donors/acceptors in region specified by "All"

(which means that in all the system) within the distance cutoff = 3 Å and angle cutoff = 135°".

The results are saved in three output files: (i) "All.UU.avg.dat" file for solute-solute H-bond averages, (ii) "All.UV.avg.dat" file for solute-solvent H-bond averages and (iii) "All.bridge.avg.dat" file for solvent bridging. It should be noted that all solvent residues are named "WAT" as the Amber default.

Furthermore, we can get a representative snapshot with VMD of a ligand showing the H-bond interactions with the surroundings (solvent molecules and protein residues).

A-2.5 Extracting snapshots with Cpptraj

After ensuring that the classical MD simulation is carried out properly, the next step is to extract some (one hundred) snapshots along the MD simulation (stabilised region in Figure 70) and take them as a starting geometries for QM/MM calculations to compute the absorption and emission spectra. The extraction of a statistical number of snapshots can be achieved using With these commands we extract snapshots from the frame range of 500 to 1000. We extract 1 snapshot each 5 frames, given 100 snapshots. It is not recommended to extract snapshot from the equilibration region (Figure 70). Now, we have one hundred snapshots stored in pdb format. It remains to convert these "snapshot.pdb" files to "Tinker.xyz" files as the QM/MM calculation uses a coupling scheme 139 between Gaussian09 138 and Tinker 140 to describe the entire system (the ligand is described at the QM level and the rest is treated at the MM level within a not polarisable model).

For each converted snapshot, solvent molecules within a radius of 15 Å from the ligand (QM subsystem) are selected and only the QM subsystem and the surrounding residues within a radius of 9 Å are allowed to move (Figure 73), given new Tinker files "Tinker.xyz, Tinker.pdb and Tinker.key". Thus, we will have three Tinker files for each snapshot for performing QM/MM calculation.

The "Tinker.xyz" file contains the Cartesian coordinates, atom types and connectivity of the system. While, the "Tinker.key" file contains the solvent molecules and the surrounding residues within a radius of 9 Å from QM subsystem. In addition, it contains the force field parameters that should be used during calculation. Thanks to the micro-iterations technique, 141 we converge the MM subsystem for every optimisation step of the QM subsystem. Once the optimisation calculations are achieved, we proceed to the calculations of the vertical transition energies using TDDFT method. Finally, for simulating the spectrum (Figure 75) we perform a Gaussian function convolution of the vertical transition energies of 100 snapshots, using a HWHM of 0.2 eV (adjustable depending of suitable spectral shape you are looking for). HWHM indicates the width of peaks for the spectrum and therefore increasing its value provides wider peaks.

Figure 1

 1 Figure 1 Scheme of firefly's bioluminescent reaction. The chemical reaction takes place inside the luciferase cavity.

P a g e | 13 Figure 2

 132 Figure 2 The six possible chemical forms of firefly oxyluciferin resulted from phenol/enol deprotonation and keto/enol tautomerization.

Figure 3

 3 Figure 3 Chemical structures of firefly oxyluciferin and their corresponding synthesised analogues (Figure reproduced from ref. 20 ). The principal moieties of oxyluciferin structure is given. The phenol-cycle presents

Figure 4

 4 Figure 4 Types of luminescence and their excitation modes (source of energy).

Figure 5

 5 Figure 5 Jablonski diagram illustrating various processes leading to luminescence. The VR indicates the vibrational relaxation, IC the internal conversion, ISC the intersystem crossing, and NRT the non-radiative transition. The lowest horizontal black line represents the ground vibrational state of the ground electronic state S0. The upper horizontal black lines represent the ground vibrational state of the excited electronic states (S1, S2, etc.). The blue thinner horizontal lines represent the vibrational levels of the electronic states.

Figure 6

 6 Figure 6 Marine and terrestrial bioluminescent species.

P a g e | 22

 22

Figure 7

 7 Figure 7 General mechanism of bioluminescence.

P a g e | 23 Figure 8

 238 Figure 8 Applications of Bioluminescence. Lighting system for urban purpose (a, b and e). Toxicants detecting test (c). Detection of infectious diseases in living mice by Bioluminescence imaging (d).

Figure 10

 10 Figure 10 General formula of amino acid.

Figure 11 a

 11 Figure 11 a) General structure of polypeptide. The ends of the peptide chain with the -COOH and -NH2 groups are known as the C-terminal and the N-terminal, respectively. b) Ribbon representation of the North American firefly "Photinus pyralis" luciferase. The small C-terminal and the large N-terminal are shown in grey and green colours, respectively. The active site of the luciferase is presented in red.

  factor. In this determinant, the N electrons are occupying N spin orbitals (𝜒 1 , 𝜒 2 , … , 𝜒 𝑁 ) without specifying which electron is in which orbital. It should be note that, the rows and the columns of an N-electron Slater determinant are labelled by electrons (first row (1), second row (2), etc.) and spin orbitals (first column χ1, second column χ2, etc.), respectively. Interchanging the coordinates of two electrons is equivalent to interchanging two rows in the determinant, which changes its sign. Thus, the Slater determinant satisfies the requirement of the anti-symmetry principle. Moreover, possessing two electrons occupying the same spin orbital of the determinant equivalent to possess two columns of the P a g e | 31

Figure 12

 12 Figure 12 Schematic representation of the molecular orbitals (MOs) classification in CASSCF calculation:

𝑡) Eq. 57 1

 571 This equation defines the TD exchange-correlation potential (𝑣 𝑥𝑐 [𝜌](𝑟, 𝑡)). The first term (𝑣 𝑒𝑥𝑡 (𝑟, 𝑡)) is the TD external potential. The second term is TD Hartree potential that accounts for the classical electrostatic interaction between electrons. The exchange-correlation P a g e | 45 and the Hartree potentials are functionals of the density, so the TD-KS equation has to be solved self-consistently. The exact expression of exchange-correlation potential (𝑣 𝑥𝑐 ) as a functional of the density is unknown. The quality of the results will obviously depend on the quality of the applied approximation to the 𝑥𝑐 potential. Thanks to the TD-DFT development, it has been possible to perform quantum mechanics calculations in the excited state of relatively large systems and obtain their excited state energies. Molecular mechanics (MM) Very large biological systems like proteins are studied mainly with molecular mechanics (MM) methods, because other methods such as quantum mechanical methods, based on approximate solution of the Schrödinger equation: semi-empirical, ab initio and DFT, are not feasible for large systems, and they are time-consuming. Molecular mechanics (MM) methods, also known as force field methods, are widely applied for calculating geometries and energies of large molecules such as proteins, peptides, nucleic acids, etc. However, MM method ignores motion of the electrons and hence, quantum effects are not considered in calculations. MM method is based on a mathematical model of a molecule as a collection of charged balls (representing the atoms) held together by springs (representing the chemical bonds).

Figure 13

 13 Figure 13 Schematic representation of the various terms that contribute to a typical force field.

  solvent's experimentally obtained dielectric constant, which generates a reaction field in response to the solute electron density. One of the first implicit solvent model, known as Born model,[START_REF] Born | [END_REF] describes a solute placed in spherical cavity surrounded by the dielectric continuum where the solute and the solvent would interact based upon the net charge. The problem with the Born model is that it did not allow for mutual polarisation between the solute and the solvent. The Born model was refined by the Onsager model 64 in which the solvent response is computed by treating the solute as a simple dipole moment in a spherical cavity surrounded by continuum dielectric. The Onsager model is appropriate for many more species. In this model, the solute dipole induces a dipole in the surrounding medium which in turn generates a reaction field in the cavity. The reaction field then interacts with the solute dipole, giving additional stabilisation for the system. The reaction field approach can be incorporated into quantum mechanical calculation where it is widely known as the self-consistent reaction field (SCRF) model. More recent models place the solute into a shaped cavity constructed from a series of spheres centred on each atom, the solute charge is then represented as a series of charge points spread out on the resulting cavity surface. The widely used model to deal with implicit solvation is the Polarisable Continuum Model 65 (PCM). The PCM model is based upon the idea of generating multiple overlapping spheres around each of atom within the solute inside of a P a g e | 50 dielectric continuum. In the PCM model, the solvation free energy (𝛥𝐺 𝑠𝑜𝑙 ) of a solute molecule can be expressed as the sum of four terms, 𝛥𝐺 𝑠𝑜𝑙 = 𝛥𝐺 𝑒𝑙𝑒𝑐 + 𝛥𝐺 𝑟𝑒𝑝 + 𝛥𝐺 𝑑𝑖𝑠 + 𝛥𝐺 𝑐𝑎𝑣 Eq. 65where the first three terms (𝛥𝐺 𝑒𝑙𝑒𝑐 , 𝛥𝐺 𝑟𝑒𝑝 𝑎𝑛𝑑 𝛥𝐺 𝑑𝑖𝑠 ) represent the electronic, Pauli repulsion, and dispersion interactions, respectively, between the solute and the implicit solvent molecules.The last term (𝛥𝐺 𝑐𝑎𝑣 ) represents the free energy of cavitation required to create a cavity of the molecular shape of the solute.The cavity employed in the PCM model is generated by a series of overlapping spheres defined by van der Waals radii of the individual atoms. There are two surfaces, solvent accessible surface (SAS) and solvent excluded surface (SES). In the SAS, the surface corresponds to the point of contact between the spheres of the solute and the sphere of the solvent. SES is the topological boundary of the union of all possible solvent molecules that do not overlap with the solute.

Figure 14

 14 Figure 14 Illustration of implicit and explicit models. On the left: PCM cavity of a solute molecule, the polarizability is represented by the bulk dielectric constant (𝜀). Implicit model accounts mainly for electrostatic interactions. On the right: A solute molecule surrounded by explicit water molecules. Explicit model accounts for specific interactions such as hydrogen bonding.

( 67 Figure 15

 6715 Figure 15 Illustration of the non-equilibrium and equilibrium solvation states. The geometry optimisation and the excitation (absorption) steps are shown.

Figure 16

 16 Figure 16 Schematic representation of a cubic periodic boundary conditions in a 3-Dimentional system. An initial molecule (in red) is represented in the central cubic simulation box (in grey), which in turn is replicated throughout the space to form an infinity lattice (white boxes).

Figure 17

 17 Figure 17 Schematic illustration of cubic periodic boundary condition in a 2-Dimentional system showing a spherical cut-off (in green line). When calculating the forces acting on a particle of interest (in red), only the particles inside the spherical cut-off are considered.

  as hydrogen bonds through the attraction between the electronegative oxygen atom and the electropositive hydrogen atoms. Several water models have been developed in order to find an accurate description of the charge distribution of the water molecules that can adequately consider the hydrogen bonding in solution. The most commonly used in MD simulations are the rigid non-polarisation models in which water molecules are represented as a set of point charges at fixed positions relative to the oxygen nucleus. The potential of these water models is composed of Lennard-Jones (for dispersion and repulsion forces) and Coulombic (for the electrostatic interactions) terms. These rigid non-polarisation models are also known as the transferable intermolecular potential with X points (TIPXP) models (e.g., TIP3P, TIP4P, SPC and SPC/E models). The TIP3P model is considered the simplest water model as it has three interaction sites corresponding to the three atoms of the water molecule. In this model, each atom is assigned a point charge and oxygen atom gets Lennard-Jones parameters. The TIP3P model is widely used in biomolecular systems simulations due to its implementation simplicity and low computational cost, and it is the one used in this Ph.D. Thesis. Furthermore, TIP3P model has been implemented in the AMBER program and most biomolecular MM force fields have been developed to be compatible with this model. P a g e | 55

Figure 19

 19 Figure19 Illustration of minimisation energy using the steepest descent algorithm. Here, the minimum is reached from the starting point (𝒙 0 ) after three (𝒙 3 ) optimisation steps, with the steps of the steepest descent method in green.

  respectively. 𝑟 𝑖 (𝑡 + 𝛿𝑡), 𝑎 𝑖 (𝑡 + 𝛿𝑡) and 𝑣 𝑖 (𝑡 + 𝛿𝑡) are the position, acceleration and velocity of atom 𝑖 at time 𝑡 + 𝛿𝑡, respectively.Due to the numerical stability, reasonable accurate and simplicity these algorithms are the most preferred choice as integrators for an MD simulation. The choice of a correct time step is important. A too large time step will cause atoms to move too far along a given trajectory (will cause large fluctuations), thus poorly simulating the motion. A too small time step will improve the accuracy of numerical solution of the classical equation of motion, but will make it necessary to run more iterations, thus taking longer to run the MD simulation. In any case time step should not be larger than the timescale of the shortest motion (intramolecular vibrational motions or time between collisions). This can give a time step on the order of five femtoseconds for simulating a Lennard-Jones fluid or rigid molecules and one femtosecond for simulating vibrating molecules.In this Ph.D. Thesis, the molecular dynamics (MD) simulations have been performed with a time step of 2 fs. The snapshots (which represent the system at each time 𝑡) obtained from MD trajectory are used for the computation of vertical excitation energies and simulation of the absorption and emission fluorescence spectra. P a g e | 61

Figure 22

 22 Figure 22 Representative pattern of the different steps of the umbrella sampling. The reaction coordinate is separated on several windows. The obtained information is then pieced together to produce an unbiased histogram. The final PMF profile is computed using WHAM program. The protein is represented in ribbon and the ligand in red inside of cavity of the protein.

Figure 23

 23 Figure 23 Docking of small molecule ligand with a protein receptor to produce a complex.

52

 52

P a g e | 65

 65 𝐸 𝑡𝑜𝑡𝑎𝑙 𝑄𝑀/𝑀𝑀 = 𝐸 𝑀𝑀+𝑄𝑀 𝑀𝑀 𝑙𝑒𝑣𝑒𝑙 + 𝐸 𝑄𝑀 𝑄𝑀 𝑙𝑒𝑣𝑒𝑙 -𝐸 𝑄𝑀 𝑀𝑀 𝑙𝑒𝑣𝑒𝑙 Eq. 76 where the subscripts QM and MM represent the atoms in the QM and MM subsystems, respectively, and the superscripts indicate the levels of theory at which their energies are computed. The commonly used subtractive QM/MM coupling scheme is the ONIOM method (Our own N-layered Integrated Molecular Orbital and Molecular Mechanics) implemented in Gaussian program and is illustrated in Figure 25.

Figure 25

 25 Figure 25 Illustration of the subtractive QM/MM coupling scheme.

Figure 26

 26 Figure 26 Illustration of the subtractive QM/MM coupling scheme.

  which the polarisation effects are included. In an electrostatic embedding model, the QM-MM electrostatic interactions between the subsystems are treated during the calculations of the electronic wavefunction. The MM point charges are incorporated as one-electron operators in the QM Hamiltonian, which is thus contains an additional operator (in atomic units): 105 𝐻 ̂𝑒𝑙𝑒𝑐 𝑄𝑀-𝑀𝑀 = ∑ ∑ 𝑍 𝐴 𝑞 𝐵 |𝑅 𝐴 -𝑅 𝐵 | 𝐵𝜖𝑀𝑀 𝐴 -∑ ∑ 𝑞 𝐵 |𝑟 𝑖 -𝑅 𝐵 | 𝐵𝜖𝑀𝑀 𝑁 𝑖 Eq. 78 P a g e | 67

Figure 27

 27 Figure 27 The link atom approach for the boundary treatment. The Q1 and M1 represent the first atom of the QM and MM subsystems, respectively. The link atom (in red) is positioned at an appropriate distance along the QM/MM bond vector.

Figure 28

 28 Figure 28 Illustration of localised orbitals method, local self-consistent field (LSCF) method, for the boundary treatment.

Figure 29

 29 Figure 29 UV-VIS absorption spectra of firefly oxyluciferin obtained by computing excitation energies in water. a) UV-VIS stick spectrum and b) convoluted UV-VIS spectrum by Gaussian function.

Figure 30

 30 Figure 30 Schematisation of the Franck-Condon (FC) principle and Herzberg-Teller (HT) approximation. FC principle is represented by the vertical dashed line showing the highest probability of transition. HT approximation represented by the grey rectangle showing the broader possibilities of transitions. The GVS indicates the ground vibrational state, EVS the excited vibrational state.

P a g e | 74 relative

 74 line shapes. VG model assumes that both initial and final PESs are equal and computes only the force constant of the initial state. It ignores the normal modes and frequencies from the final state and only evaluates the energy gradient of the final state about the equilibrium structure of the initial state.

Figure 31

 31 Figure 31 Vibrationally-resolved absorption spectra of firefly natural oxyluciferin generated by using FC, HT and FCHT approximations in water (PCM) within TI framework and AH model. The emitted intensity by second is given in µJ/mol, before being normalised to one.

Figure 32

 32 Figure 32 Representation of the initial (in blue) and final (in red) state geometries of the natural oxyluciferin (left) and its analogue (right) obtained after the superposition procedure within AH model. The initial and final state geometries of Natural oxyluciferin (left) is superposed correctly. The initial and final state geometries of the oxyluciferin analogue (right) is superimposed, which explains why the spectrum cannot be generated.

Figure 33

 33 Figure 33 Vibrationally-resolved absorption spectra computed with the FCHT approximation in gas phase for Phenol-cycle analogue (a) protonated, b) deprotonated) using the Adiabatic Hessian (AH) model. In the case of the protonated phenol-cycle analogue, the superimposition of the electronic ground and first excited states minimum energy structures is provided and explains why the spectrum cannot be correctly simulated (The spectrum calculation has been forced).

Figure 34

 34 Figure 34 Vibrationally-resolved absorption spectra computed with the FCHT approximation in water (PCM) for a) protonated phenol-OMe (left: with VG model, right: with AH model) and b) protonated OMeenol (left: with VG model, right: with AH model).

FCHT

  within AH model. We found that the spectral shapes obtained with VG and AH models for the description of the initial and final PESs are very different. The spectra simulated with AH model represent a larger and more intense amount of vibronic transitions in the spectral region compared to the VG spectra. Furthermore, in both AH and VG models, the predicted 0-0' vibronic transition corresponds to the most intense one. Moreover, the 0-0' bands predicted with VG model are blue-shifted around 0.05 eV comparing with those obtained with AH model.The VG model is computationally cheaper than AH model, but it neglects the difference in the Hessians of the initial and final states (VG computes only the force constant of the initial state) and only evaluates the energy gradient of the final state about the equilibrium structure of the initial state, which can impair the reproduction of the spectrum. However, VG model offers the possibility to generate vibronic spectrum for a more flexible macromolecular system, which is not infeasible with AH model. In this Ph.D. Thesis, we have used both FC and FCHT approximations applying both AH and VG models within time independent framework.

Figure 35

 35 Figure 35 Simulated emission spectra updating the parameters after the first, second and third iteration for the oxyluciferin phenolate-cycle analogue. (Figure reproduced from ref. 20 ).

Figure 36 a

 36 Figure 36 a) Simulated emission spectra of the oxyluciferin phenolate-keto form considering an increasing number of MD snapshots. b) Simulated emission spectra of oxyluciferin phenolate-cycle analogue with optimising or not its geometry in the excited state. (Figure reproduced from ref. 20 ).

Firefly's luciferase

  contains two terminal domains: N-terminal and C-terminal containing 400-500 and 110-130 residues, respectively. The luciferase active site is located at the interface between the two terminals. During the bioluminescence process the luciferase can adopt two different conformations by alternating the C-terminal. The first luciferase catalytic conformation is characterised by an open form catalysing the initial adenylation reaction of the luciferin (the step 1 in Figure 1). The second catalytic conformation adopts a closed form caused by a rotation of the C-terminal (about 140°) during the last reaction leading the light emission. P a g e | 82

Figure 38

 38 Figure 38 Graphical representation of a) 5'-O-[N-dehydroluciferyl-sulfamoyl]-adenosine (DLSA), b) oxyluciferin phenolate-keto and c) mono-deprotonated adenosine monophosphate (AMPH) structures.

  A representative MD snapshots of phenolate-keto, phenolate-enol and phenol-enolate forms showing the H-bonding between the heteroatoms and water molecules are presented in Figure 39. For instance in the representative MD snapshots of phenolate-keto, phenolate-enol and phenol-enolate (Figure P a g e | 88

Figure 40

 40 Figure 40 Normalised simulated absorption spectra (continuous lines) of both the protonated (red lines) and deprotonated (blue lines) forms and the normalised experimental absorption spectra (dashed lines) recorded at pH 5 (red lines) and pH 11 (blue lines) of A) phenol-cycle, B) phenol-OMe, and C) OMe-enol. (Figure reproduced from ref. 20 ).

Figure 41

 41 Figure 41 Comparison of the simulated absorption spectra of the each analogue (red lines) and its corresponding oxyluciferin form (blue lines) simulated with QM/MM methods considering 100 MD snapshots. S1, S2 and S3 are taken from the MD snapshot that corresponds to the spectrum maximum (Figure reproduced from ref. 20 ).

Figure A- 1 ,

 1 Figures A-X refer to figures in the Appendix)

Figure 42

 42 Figure 42 Simulated absorption spectra and molecular orbitals involved in each electronic transition of A) phenol-cycle, B) phenol-keto. (Figure reproduced from ref. 20 ).

(

  deprotonated). The results show that the band characterising the emission of the protonated form rapidly decays with the rise of the one characterising the deprotonated form indicating that the excited state dynamics of the phenol-cycle analogue is ruled by an ultrafast ESPT reaction.

Figure 43

 43 Figure 43 Normalised simulated emission spectra (continuous lines) of both the protonated (red lines) and deprotonated (blue lines) forms and the normalised experimental emission spectra (dashed lines) recorded at pH 5 (red lines) and pH 11 (blue lines) of A) phenol-cycle, B) phenol-OMe, and C) OMe-enol. (Figure reproduced from ref.20).

Figure 44

 44 Figure 44 Comparison of the simulated emission spectra of each analogue (blue lines) and its corresponding oxyluciferin form (green lines) simulated with QM/MM methods considering 100 MD snapshots. (Figure reproduced from ref. 20 ).

P a g e | 95 P a g e | 96 P a g e | 97 P a g e | 98 P a g e | 99 P

 9596979899

Figure 45

 45 Figure 45 Simulated absorption spectra of the oxyluciferin natural forms, in the gas phase (grey), with the LR-PCM formalism (red) and with QM/MM methods considering 100 MD snapshots (blue). (Figure reproduced from ref. 135 ).

Figure 46

 46 Figure 46 Simulated emission spectra of the oxyluciferin natural forms, in the gas phase (grey), with the LR-PCM formalism (red) and with QM/MM methods considering 100 MD snapshots (blue). (Figure reproduced from ref. 135 ).

Figure 47

 47 Figure 47 Comparison of the FC and FCHT vibronic absorption and emission spectra of phenolate-keto computed a) in the gas phase (GP), b) with implicit model (PCM), c) with microsolvation in GP (µSolv-GP) and d) with microsolvation in PCM (µSolv-PCM). (Figure reproduced from ref. 135 ).

Figure 48 a

 48 Figure 48 a) absorption and b) emission vibronic spectra of phenolate-keto in gas phase (GP), in implicit solvent (PCM) and in microsolvated systems: in GP (microsolv-GP) and in PCM (microsolv-PCM). The vibronic transitions of the dominant-intense band (grey shadow) and the medium-intense band (beige shadow) are given as sticks (left side for absorption and right side for emission). The vibronic spectra have been shifted to match the spectra maxima obtained with the PCM model (Figure reproduced from ref. 135 ).

P a g e | 119 P

 119

Figure 49

 49 Figure 49 Synthetic oxyluciferin analogues and doubly protonated adenosine monophosphate (AMP) structures.

Figure 50

 50 Figure 50 Vibronic absorption spectra of the oxyluciferin analogue forms computed in gas phase using FCHT approximation with Adiabatic Hessian (AH) model. a) phenol-cycle, b) phenol-OMe and c) OMe-enol; left side: protonated form, right side: deprotonated form. The vibronic transitions are depicted with sticks. The calculations of the protonated phenol-cycle spectrum have been forced as the excited state structure of this compound is distorted with respect to the one of the ground state, giving unreliable spectrum. (Figure reproduced from ref. 150 ).

P a g e | 134 Figure 51

 13451 Figure 51 Experimental (Solid line) and computed (Dashed line) normalised absorption spectra of a) phenol-cycle, b) phenol-OMe and c) OMe-enol in water at pH 5, 7 and 11. The deprotonated forms populations are calculated out of the CpHMD trajectories and are indicated as "% deprotonated". The absorption wavelength maxima obtained in water and the ones considering or not the vibrational information computed in gas phase (GP) are given in nm and eV. (Figure reproduced from ref. 150 ).

P a g e | 135 Figure 52

 13552 Figure 52 Vibronic absorption spectra computed with the FCHT approximation in water (PCM) for a) Phenol-OMe analogue (left: protonated, right: deprotonated) and b) OMe-enol analogue (left: protonated, right: deprotonated) using the Adiabatic Hessian (AH) model. λmax,abs is the wavelength of absorption spectrum maximum and the values in brackets are the vertical energies computed with TDDFT using PCM.

  Figure A-2a) in order to investigate the influence of the used VG instead AH model on the vibronic calculations. In general, the spectral shapes of VG spectra are very different from the AH one. The AH spectra represent a larger and more intense amount of vibronic transitions in the spectral region around 390-450 nm (Figure 52), while in the VG spectra model these transitions (around 410-430 nm) are less intense (Figure 53a and Figure A-2a). This could explain the broader spectral shape observed using the AH model. Furthermore in both AH and VG model, the predicted 0-0' vibronic transition is the most intense one, except for the predicted one of phenol-OMe obtained with AH model but, it still closer to the wavelength maximum of vibronic spectrum. Moreover, the 0-0' band predicted with VG model are blue-shifted around 0.05 eV comparing with those obtained with AH model. This difference in energy is rather small, notwithstanding the VG model weaknesses.Furthermore, the vibronic absorption spectra simulated with or without AMP using VG model exhibit similar spectral shapes, characterised by a quite intense peak corresponding to the 0-0' vibronic transition (Figure53 and Figure A-2). The fact that the 0-0' vibronic transition is by far the most intense one explains the sharp shape of the simulated absorption spectra. In addition, all the spectra show also some intense band appears as a consequence of low intense vibrational transitions (zoomed regions in Figure53 and Figure A-2).

P a g e | 138

 138 neutral analogues and their deprotonated forms in presence or not of AMP, we find mainly the same vibronic transitions. Apart from the differences discussed above on the vibronic transitions considering AMP, we observe that the most intense band of the deprotonated forms are redshifted around 107-127 nm (compounds alone) and 102-112 nm (with mono-deprotonated AMP). These redshifts of the maximum bands are less important (40-41 nm) when the doublydeprotonated AMP is considered (Figure53and Figure A-2).

Figure 53

 53 Figure 53 The vibronic spectra of phenol-OMe (left) and its deprotonated form (right) with or without AMP (with mono-or doubly-deprotonated) computed with the FCHT approximation in water using PCM model within Vertical Gradient model. a) analogue alone, b) analogue-AMPH -and c) analogue-AMP 2-. The vibronic transitions are depicted with sticks and those with less intensity (area framed by a green rectangle) are zoomed.

Figure 54

 54 Figure 54 Absorption spectra of a) phenol-cycle-AMP, b) phenol-OMe-AMP and c) OMe-enol-AMP complexes in water between pH 5 and pH 11. (Top) experimental and (Bottom) computed unnormalised spectra. The deprotonated forms populations are calculated out of the CpHMD trajectories and are indicated as "% deprotonated". The absorption wavelength maxima of the protonated (pH 5) and deprotonated (pH 9 and 10) forms are given in nm and eV. (Figure reproduced from ref. 150 ).

P a g e | 142 P

 142

Figure 55

 55 Figure 55 Natural firefly oxyluciferin structure, Natural_Oxy, and the synthesised analogues under study (structural modification depicted in blue). Their experimental bioluminescence values (in nm) are taken from ref. 35-37 (Figure reproduced from ref. 151 ).

Table 5

 5 Vertical transition energies (Te), the corresponding wavelengths (λabs and λemi) and oscillator strengths (f) of the absorption and emission of Natural_Oxy and its analogues in gas phase, in implicit water and in protein. For the compounds in protein, the Te values correspond to the wavelengths of the absorption and emission spectra maximum and the f values are taken from representative MD snapshots.

Figure A- 3 )

 3 FigureA-3). By analysing the NTOs, we found that their electron density computed in GP are quite similar to the one in PCM. In addition, we observed that the electron density of all the compound structures are transferred from the benzothiazole moiety to the thiazolone one during the absorption transition and conversely for the emission transition (Figure56and FigureA-3). If we analyse in more detail the distribution of the electron density in the Natural_Oxy and the analogue structures, we observe that it is sensitive to the structural modifications. For instance in the case of the benzothiophene_Oxy (where CH group replaces the N atom of the Natural_Oxy benzothiazole moiety), the unoccupied NTO electron density located in the CH

P a g e | 159 Figure 56

 15956 Figure 56 Natural transition orbitals (NTOs) of Natural_Oxy and its analogues computed at the excited state minima corresponding to the emission transition in PCM. (Figure reproduced from ref. 151 ).

  40° and -160°) of the dihedral angle of the allyl group during the MD simulation (Figure A-4e). The simulation of 100 vertical S1→S0 (emission) transitions energy for the Natural_Oxy and its analogues shows that the optically bright transition energies are found in a narrow energy range, except for allylbenzothiazole_Oxy analogue where these transition energies are located in a larger range of energy (Figure A-4). Hence, a Gaussian convolution of these transition energies gives a broader allylbenzothiazole_Oxy spectrum compared to the others.

Figure 57

 57 Figure 57 Emission spectra of Natural_Oxy and its synthetic analogues a) simulated with QM/MM methods and b) measured experimentally. 35-37 The experimental bioluminescence spectra are adapted with permission from (C. C. Woodroofe, P. L. Meisenheimer, D. H. Klaubert, et al., Novel Heterocyclic Analogues of Firefly Luciferin, Biochemistry, 2012, 51, 9807-9813). Copyright (2012) American Chemical Society, and from (S. Ioka, T. Saitoh, S. Iwano, et al., Synthesis of Firefly Luciferin Analogues and Evaluation of the Luminescent Properties, Chem. Eur. J., 2016, 22, 9330-9337. Copyright (2016) WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim. (Figure reproduced from ref. 151 ).

Figure 58 a

 58 Figure 58 a) Statistical average number of H-bonding involving specific atoms of the compounds under study, water molecules, AMPH and protein active site along MD simulation in the excited state (ES). b) Atom numbering. c) Representative MD snapshots of Natural_Oxy (and its analogues) inside protein active site. (Figure reproduced from ref. 151 ).

P a g e | 164 Figure 59

 16459 Figure 59 a) Distance (in Å) between the oxygen atoms (O10 and O11) of the compound under study and carbon atom (CX) of amino acids (or the phosphate (P) atom of AMPH) calculated from MD simulations in the ES b) Atom numbering and the initial distances (in Å) between (O10 -CX) and (O11 -P) before the Natural_Oxy MD simulation. c) The calculated distance (in Å) between the atoms O10 and CX of ARG337 amino acid from two benzothiophene_Oxy MD simulations. (Figure reproduced from ref. 151 ).

Figure 60

 60 Figure 60 NTOs of Red_1_Oxy and Red_2_Oxy conformers computed at the equilibrium S1 geometry corresponding to S1-> S0 (emission) electronic transition in gas phase (GP) and in water (PCM). (Figure reproduced from ref. 151 ).

  large dihedral angle fluctuation of the allyl group during the MD simulations as it has been found for the allyl group of allylbenzothiazole_Oxy (Figure A-4e).

Figure 61

 61 Figure 61 Emission spectra of the Red_2_Oxy conformer, Natural_Oxy and allylbenzothiazole_Oxy analogue simulated in protein with QM/MM methods considering 100 statistical MD snapshots (HWHM of 0.20 eV).

P a g e | 169 P a g e | 170 P a g e | 171 P a g e | 172 P a g e | 173 P a g e | 174 P a g e | 175 P a g e | 176 P a g e | 177 P a g e | 178 P a g e | 179 Conclusion

 169170171172173174175176177178179 During this Ph.D. Thesis, theoretical studies focused on the last step of the firefly bioluminescent reaction, that is, the electronic transition that yields the light emission, have been performed. In particular, firefly's light emitter structure, oxyluciferin, and some of its synthetic analogues have been studied in water solution to get insight into the emission mechanism and identify the chemical form of the light emitter. To clarify and get insight into the firefly's bioluminescence colour tuning, an investigation study on the influence of structural modifications within the light emitter structure on the emission spectrum has been also performed considering explicitly the water molecules and the protein environment.

Figure A- 2

 2 Figure A-2 The vibronic spectra of OMe-enol (left) and its deprotonated form (right) with or without AMP (with mono-or doubly-deprotonated) computed with the FCHT approximation in water using PCM model within Vertical Gradient model. a) analogue alone, b) analogue-AMPH -and c) analogue-AMP 2-. The vibronic transitions are depicted with sticks and those with less intensity (area framed by a green rectangle) are zoomed.

(

  Figure reproduced from ref. 151 ).

(

  S0→S1), 𝐶𝑇 = 𝑞 S1 𝑡ℎ𝑖𝑎𝑧 -𝑞 S0 𝑡ℎ𝑖𝑎𝑧 with the charges taken at the equilibrium S0 structure in PCM. For emission (S1→S0), 𝐶𝑇 = 𝑞 S0 𝑡ℎ𝑖𝑎𝑧 -𝑞 S1 𝑡ℎ𝑖𝑎𝑧 with the charges taken at the equilibrium S1 structure in PCM. Negative CT values indicates a fraction of the electronic charge transferred from benzothiazole moiety to thiazolone moiety during the transition. (Table reproduced from ref. 151 ).

Figure A- 3

 3 Figure A-3 Natural transition orbital (NTO) of Natural_Oxy and its analogues computed at the S0 and S1 equilibrium geometries corresponding to S0 -> S1 (Top: absorption) and S1-> S0 (Bottom emission) electronic vertical transitions, respectively, in gas phase, in water (PCM) and in protein. For the chemical forms in protein, the selected MD snapshots correspond to the maximum absorption wavelength. The NTOs are plotted with Avogadro package. (Figure reproduced from ref. 151 ).

Figure A- 4

 4 Figure A-4 Oscillator strengths against transition energies for the 100 vertical S1→S0 (emission) transitions of Natural_Oxy and its synthetic analogues simulated with QM/MM methods used to simulate the emission spectra. a) benzothiophene_Oxy, b) dihydropyrrolone_Oxy, c) Natural_Oxy, and d) allylbenzothiazole_Oxy. Under figure d, three geometries of allylbenzothiazole_Oxy along the MD simulation with different allyl group rotation. e) Dihedral angle fluctuation obtained with CPPTRAJ program defined by the four atoms C6, C7, C8 and C9 of allylbenzothiazole_Oxy along the MD simulation in the ES. (Figure reproduced from ref. 151 ).

Figure A- 5

 5 Figure A-5 Absorption spectra of the Natural_Oxy and its analogues simulated in protein using QM/MM methods (considering 100 statistical MD snapshots, HWHM of 0.20 eV). (Figure reproduced from ref. 151 ).

Figure A- 7

 7 Figure A-7 Statistics of water molecules found around 5 Å of the carbon atom C7 (see atom numbers inFigure A-4e) of Natural_Oxy (in green) and allylbenzothiazole_Oxy analogue (in red) during the MD simulations. (Figure reproduced from ref. 151 ).

Figure A- 6 H

 6 Figure A-6 H-bonding networks formed during MD simulations at the excited state S1 of Natural_Oxy (in green) with a) benzothiophene_Oxy (in blue), b) with dihydropyrrolone_Oxy (in light-blue) and c) with allylbenzothiazole_Oxy (in red). The extracted MD snapshots correspond to the maximum emission wavelength. (Figure reproduced from ref. 151 ).

P a g e | 190 Figure 62

 19062 Figure 62 a) Optimised geometry coordinates from a Gaussian output in gas phase calculation at the QM level. b) Geometry coordinates saved in file.xyz (21 corresponds to the number of atom of the ligand).

Figure 63

 63 Figure 63 Pdb format of the optimised geometry of a chromophore a) generated with VMD and b) modified manually (changing X by LIG and the atom names O, C by O1, C1, etc.). The second column(1, 2, etc.) represents the order number of each atom, third column (O1, C1, etc.) the atom name of the force field, the fourth column the name of the ligand (LIG) which has a classification number(1). The rest is xyz coordinates.

Figure 64

 64 Figure 64 Example of a generated "file.prep" file with Antechamber.

P a g e | 192 Figure 65

 19265 Figure 65 Force field parameters for the Ligand generated by parmchk2 program and edited manually by using gaff.dat parameters file as reference.

Figure 66

 66 Figure 66 Example of leaprc file.

Figure 68

 68 Figure 68 Input file of the heating step.

Finally

  , under NPT (isothermal-isobaric ensemble) conditions within T = 300 K and P = 1 atm, we perform a classical MD simulation of 20 ns with periodic boundary conditions using the following commands: → mpirun -np 4 sander.MPI -O -i file_MD.in -o file_MD.out -p file_box.top -r file_MD.rst -ref file_heat.rst -x file_MD.mdcrd -e file_MD.ene -c file_heat.rst & MD simulation can be also performed using another Sander version called PMEMD (stands for Particle Mesh Ewald Molecular Dynamics) module. PMEMD module offers parallel CPU and it is not time-consuming compared to the Sander module. PMEMD module is available on our computer cluster and it can be used as follows: → $LAUNCHER $LOPTS $NOPT $NCPU pmemd.MPI -O -i $WDIR/file_MD.in -o $WDIR/ file_MD.out -p file_box.top -c file_heat.rst -r file_MD.rst -ref file_heat.rst -x file_MD.mdcrd -e file_MD.ene

Figure 69

 69 Figure69 Input file of the molecular dynamics simulation step. During this MD simulation, the pressure and temperature are maintained by using the Berendsen algorithm.[START_REF] Berendsen | Molecular dynamics with coupling to an external bath[END_REF]

Figure 70

 70 Figure 70 Root mean square deviation (RMSD) analysis obtained from a classical MD simulation. The calculated mass-weighted RMSD are obtained using only the heavy atoms in the ligand, the protein residues and AMPH. The structure of the first trajectory frame (first snapshot) is taken as reference.

Furthermore, in order

  to facilitate the visualisation of the dynamic with VMD we can center the coordinates of the molecule of interest by using Cpptraj program as follows: → cpptraj → parm file_box.top → trajin file_MD.mdcrd

Figure 71

 71 Figure 71 Output file containing solute-solute hydrogen bond averages.

Figure 72 A

 72 Figure 72 A representative snapshot of a ligand showing the H-bond interactions with its surrounding environment.

Cpptraj program as follows: → cpptraj P a g e | 198 → 5 →

 1985 parm file_box.top → trajin center_file_MD.mdcrd 500 1000 trajout snapshot pdb multi → run

Figure 73

 73 Figure 73 Example of an extracted snapshot from MD simulation. The circled parts are considered in the QM/MM calculation, that is, the protein, ligand and the surrounding water molecules within a radius of 15 Å.

Figure 74

 74 Figure 74 Example of QM/MM optimisation input file. To compute the vertical transition energy we put TD (with excited state of interest) instead OPT. The last line of this input considers the Tinker files (.xyz and .key files) and mentions the number of atoms in the Tinker.xyz file.

Figure 75

 75 Figure 75 Vertical transition energies of 100 snapshots computed at QM/MM level (left side) and the spectrum obtained after Gaussian convolution of the vertical transition energies of 100 snapshots (right side).

  position vectors 𝑅 𝐴 and 𝑟 𝑖 , respectively, is given in atomic units byThe two first terms (𝑇 ̂𝑒 and 𝑇 ̂𝑁) are the kinetic energy operators of electrons and nucleus, respectively. The third term (𝑉 ̂𝑒𝑁 ) represents the coulomb attraction between an electron 𝑖 and a nucleus 𝐴. The fourth and the fifth terms (𝑉 ̂𝑒𝑒 and 𝑉 ̂𝑁𝑁 ) represent the repulsion between pairs of electron and nuclei, respectively. 𝛻 𝑖2 and 𝛻 𝐴 2 are Laplacian operators, 𝑍 𝐴 and 𝑀 𝐴 are the atomic number and mass of nucleus 𝐴, respectively. The distance between two electrons 𝑖 and 𝑗, between electron 𝑖 and nucleus 𝐴, and between two nuclei 𝐴 and 𝐵 are represented by 𝑟 𝑖𝑗 , 𝑟 𝑖𝐴 and 𝑅 𝐴𝐵 , respectively.

														P a g e | 29
				𝐻 ̂ = 𝑇 ̂𝑒 + 𝑇 ̂𝑁 + 𝑉 ̂𝑒𝑁 + 𝑉 ̂𝑒𝑒 + 𝑉 ̂𝑁𝑁			Eq. 7
	𝐻 ̂ = -	1 2	∑ 𝛻 𝑖 2 𝑛 𝑖	-	1 2	∑ 𝑀 𝐴	𝛻 𝐴 2 𝑀 𝐴	-∑ ∑ 𝑀 𝐴 𝑛 𝑖	𝑍 𝐴 𝑟 𝑖𝐴	+ ∑ ∑ 𝑛 𝑗>1 𝑛 𝑖	1 𝑟 𝑖𝑗	+ ∑ ∑ 𝑀 𝑀 𝐵>𝐴 𝐴	𝑍 𝐴 𝑍 𝐵 𝑅 𝐴𝐵	Eq. 8

.6.1 Basis sets classification

  𝑔(𝛼, 𝑙 𝑥 , 𝑙 𝑦 , 𝑙 𝑧 , 𝑥, 𝑦, 𝑧) = 𝑁𝑒 -𝛼𝑟 2 𝑥 𝑙 𝑥 𝑦 𝑙 𝑦 𝑧 𝑙 𝑧 Eq. 33

	P a g e | 35
	type orbitals (GTOs) are obtained from a radial part and a spherical part. The GTOs have a form
	depending on Cartesian coordinates:
	Here, the radial part has a different dependence on 𝑟. The exponents 𝑙 𝑥 , 𝑙 𝑦 , 𝑙 𝑧 in Cartesian
	coordinates should not be confused with the quantum numbers. The sum of these exponents
	(𝐿 = 𝑙 𝑥 + 𝑙 𝑦 + 𝑙 𝑧 ) is used to define the angular momentum of basis functions. The Gaussian type
	function (GTF) is known as: s-type (𝐿 = 0), p-type (𝐿 = 1) and d-type (𝐿 = 2), etc.
	2.1Basis set can be classified into the following types:
	Minimal basis sets, such as STO-3G, STO-4G and STO-3G*, contain minimum number of
	basis functions that are needed for each atom. The STO-3G, STO-4G and STO-3G* basis sets use
	three and four gaussian primitives (3G and 4G) per basis function. The STO-3G* is the polarised
	(with signe *) version of STO-3G.
	Pople basis sets, such as 3-21G, 6-31G, 6-311G and 6-311+g* (with polarised* and
	diffuse + functions), contain more basis functions per atom than the previous one. This type of
	basis sets have two or more sizes of basis function for each valence orbital. For instance, the 6-
	31G basis set is a double split valence basis set in which the core orbital is described by one of
	contracted Gaussian function comprised of 6 primitives Gaussians. The valence orbital is split
	into two sets of contracted Gaussian functions, 3 and 1 primitive Gaussians. The 6-311G basis
	set is a triple split valence basis set in which the core orbital is a contraction of 6 primitive
	The Gaussian function

  Here 𝐹 [𝜌] is a universal HK functional. It operates only on density 𝜌 and it is independent from the external potential 𝑣(𝑟). 𝑉 𝑛𝑒 is the potential energy of nucleus-electron attraction, 𝑇 is the electron kinetic energy and 𝑉 𝑒𝑒 is the electron-electron repulsion energy.

							P a g e | 41
	𝑉 𝑛𝑒 = ∫ 𝜌(𝑟) 𝑣(𝑟)𝑑𝑟	;	𝐹 [𝜌] = 𝑇 [𝜌] + 𝑉 𝑒𝑒 [𝜌]	Eq. 42
	The second HK theorem pointed out that there is an equivalent to the variational
	principle in molecular orbital theory, that is, any trial ground state electron density 𝜌 ′ will give
	an energy functional 𝐸[𝜌′] that is greater or equal than the real ground state energy 𝐸 0 of the
	system with density 𝜌 0 .
			𝐸[𝜌′] ≥ 𝐸 0 [𝜌 0 ]	Eq. 43
	In order to obtain the best density, the energy has to be minimised with respect to
	density variations, just as one does in determining the best wavefunction 𝜓 in variational
	molecular orbital theory. The total ground state energy will reach the minimum uniquely if the
	electron density corresponds to the exact ground state electron density.
	Another approach is developed to calculate the energy directly without recourse of the
	wavefunction, named as Kohn-Sham (KS) approach. The premise of the KS approach is to
	introduce an auxiliary system of non-interacting electrons into the formalism. This leads to the
	central equation in Kohn-Sham DFT, which has the form of a single particle Schrödinger
	equation:
	[-	1 2	𝛻 𝑖 2 + 𝑣(𝑟) + ∫	𝜌(𝑟 ′ ) |𝑟 -𝑟 ′ |	𝑑𝑟 ′ + 𝑣 𝑥𝑐 (𝑟)] 𝜙 𝑖 = 𝜀𝜙 𝑖	Eq. 44
	functional of the electron density 𝜌(𝑟):
	𝐸[𝜌] = 𝑉 𝑛𝑒 [𝜌] + 𝑇 [𝜌] + 𝑉 𝑒𝑒 [𝜌] = ∫ 𝜌(𝑟) 𝑣(𝑟)𝑑𝑟 + 𝐹 [𝜌]	Eq. 41
	where

  In this case, the velocity is explicitly computed at half-time steps, which in turn is used to update the position at full time step.

							P a g e | 60
	𝑣 𝑖 (𝑡 +	1 2	𝛿𝑡) = 𝑣 𝑖 (𝑡 -	1 2	𝛿𝑡) + 𝑎 𝑖 (𝑡)𝛿𝑡	Eq. 72
	𝑟 𝑖 (𝑡 + 𝛿𝑡) = 𝑟 𝑖 (𝑡) + 𝑣 𝑖 (𝑡 +	1 2	𝛿𝑡) 𝛿𝑡	Eq. 71

1 2 𝛿𝑡.

  If the QM and MM subsystems are not covalently bound to each other, one can use the previous models to describe the system by defining the subsystems such that the boundary does not pass through a covalent bond. On the other hand, if a covalent bond exists at the QM/MM boundary, special care is required when evaluating the QM wavefunction. A straightforward cut of the chemical bond connecting the QM and MM subsystems can create one or more unpaired

electrons in the QM subsystem and hence, cause unsaturation in the QM subsystem. These electrons are paired in bonding orbitals with electrons that belong to the atom on the MM P a g e | 68

Table 1

 1

), we can classify the chemical forms of the oxyluciferin and the analogues in three groups: (i) group with a low charge transfer (LCT group) where the CT character is almost negligible (|qS1 thiaz -qS0 thiaz |≤ 0.1e), (ii) group with a high charge transfer (HCT group) where the CT character is equal or larger than 0.25e and (iii) group with a medium charge transfer (MCT) presenting an intermediate CT character. Following these indications, phenolenol and its analogue can be classified to the LCT group as the CT character of the transition is quite small. Phenolate-keto, phenolate-enol, phenolate-enolate and their analogues, presenting an intermediate CT character, can belong to the MCT group. Phenol-keto and phenol-enolate forms and their corresponding analogues can be classified to the HCT group as they show a large CT character.

Table 3

 3 Stokes shifts (eV) (𝐸 𝑡ℎ𝑒𝑜 𝑎𝑏𝑠 -𝐸 𝑡ℎ𝑒𝑜

𝑒𝑚𝑖 ) for oxyluciferin forms (coloured lines) and their analogues (uncoloured lines) computed in gas phase and with different solvation models. Experimental values (𝐸 𝜆𝑚𝑎𝑥,𝑒𝑥𝑝 𝑎𝑏𝑠 -𝐸 𝜆𝑚𝑎𝑥,𝑒𝑥𝑝 𝑒𝑚𝑖 ) are also given. (Table reproduced from ref. 135 ).

Table A -

 A

2 Square of the transition dipole moment (µge in a.u.) between the ground and excited states of the absorption and emission in gas phase and in water (PCM). (Table reproduced from ref. 151 ).

Table A -

 A 5 Mulliken atomic charges of benzothiazole and thiazolone moieties at ground (S0) and first (S1)excited states computed at the equilibrium S0 and S1 structures in PCM, corresponding to the absorption and emission minimum energy structures, respectively. The CT character is calculated as the difference of the thiazolone moiety charges between the final state 𝑞 Sf 𝑡ℎ𝑖𝑎𝑧 and initial state 𝑞 Si 𝑡ℎ𝑖𝑎𝑧 . For absorption

  It is also possible, with Cpptraj program, to track and find not only the solute-solute Hbonds, but also solute-solvent H-bonds by using the following commands: All out All.hbvtime.dat solventdonor :WAT solventacceptor :WAT@O avgout All.UU.avg.dat solvout All.UV.avg.dat bridgeout All.bridge.avg.dat → run

		P a g e | 197
	A-2.4.2	Solute-solvent and solute-solute H-bonds
	→ cpptraj
	→ parm oxy_box_last.top
	→ trajin md1_center.mdcrd
	→ hbond
	protein-ligand interactions.

Acknowledgments

Appendix

Table A-1 Vertical transition energies (Te), the corresponding wavelengths (λabs and λemi) and oscillator strengths (f) of the absorption and emission transitions of phenolate-keto in water (PCM) computed using different functionals (B3LYP, CAM-B3LYP, M06, M06-2X and PBE0) with 6-311g(2d,p) basis set. The computational cost to obtain the transition energies is evaluated as the CPU days needed for optimisation and transition energies calculation. The experimental results have been taken from ref. 20

A-2.2 Running classical MD simulation with Sander

Once the solvation box is correctly set and the input files file_box.top and file_box.crd for Sander module of the Amber14 program are generated, the next step is to run a classical molecular dynamics (MD) simulation with Amber14 program 123 using an Amber (16) module which performs the energy minimisation and molecular dynamics. This module is known as Sander module (acronym for Simulated Annealing with NMR-Derived Energy Restraints). First, we start by carrying out the energy minimisation by using the following commands

→ module load ambermpi → mpirun -np 4 sander.MPI -O -i file_min.in -o file_min.out -p file_box.top -r file_min.rst -ref fle_box.crd -x file_min.mdcrd -e file_min.ene -c file_box.crd &

In this minimisation step, all the system is allowed to move freely along the simulation, except for the ligand/chromophore which is restrained by applying some restraint. The input file "file_min.in" looks like: After the energy minimisation we take the "file_min.rst" as a restart point to run the heating step. In this step the system is heated from 100 K to 300 K in 20 pico-seconds while restraining the ligand:

→ mpirun -np 4 sander.MPI -O -i file_heat.in -o file_heat.out -p file_box.top -r file_heat.rst -ref file_min.rst -x file_heat.mdcrd -e file_heat.ene -c file_min.rst &

The input file "file_heat.in" looks like: (16) For more information regarding the AMBER commands, check out the following AMBER website: https://ambermd.org/Manuals.php.