Can We Describe Biological Systems with Quantum Mechanics?

Quantum Mechanics is the favourite theory to predict the structure of any group of atoms, including biological molecules. Due to numerous difficulties, however, it is necessary to introduce a series of approximations to overcome such impediments. We present a coarse-grained model of circular dichroism (CD) that is based on the theory of optical activity, developed by DeVoe, in order to predict CD spectra. In first stage, we determine the polarisability of individual monomers (residues, in the case of peptides) from experiments of molar absorptivity. The complex polarisabilities are used together with peptide structures obtained by density functional theory and other methods to determine their corresponding CD spectra, which are in reasonable agreement with their experimental counterparts.


Introduction
The description and full understanding of biological systems is one of the major scientific tasks for this century. For example, malaria killed about half a million people around the world in 2015 and that is good news because it used to kill twice as many persons, 20 years earlier. In 2015, Youyou Tu received the Nobel prize in medicine for using a plant extract, artemisinin, in the treatment of malaria. In spite of important efforts, there is still no vaccine for malaria and the complete cure remains an open problem as in the case of many other diseases. The comprehension of biological systems can be a challenge because they are composed of a variety of molecules such as lipids, proteins, and nucleic acids like DNA and RNA. These molecules have specific functions that include interaction with other molecules and, eventually, they can modify their shape to perform a certain activity.
The typical biological systems are characterised by experimentally accessible distances (from nanometers to microns) and relatively low energies (of the order of the thermal energy, k B T ). Therefore, the preferred approach to describe such systems should be Quantum Mechanics (QM). For instance, we can propose the solution of the time-independent non-relativistic Schrödinger equations for a system consisting of M nuclei and N electrons that mainly interact through a Coulomb potential. A reduced description, that captures the intended behaviour of the system, requires more than 10 6 atoms (including water molecules). Thus, it is basically impossible to solve the resulting highly-coupled Schrödinger equations, even at the level of the Born-Oppenheimer approximation. The exclusion of the solvent molecules provides a significant reduction of the number of atoms, but it is still necessary to utilise an appropriate theoretical framework. For this purpose, we briefly describe the Density Functional Theory (DFT) in section 2.
In this work, we used DFT to obtain the two stable configurations of a relatively small peptide, namely, polylysine of 150 residues. Such structures are then utilised as a reference in a proposed model of circular dichroism (CD), which is an experimental technique that takes advantage of the chiral nature of biological molecules, allowing the identification of structural changes. Once calibrated, we make use of additional approximations in order to compare the results of our model with experimental data, as explained in the following sections.

Theoretical background 2.1. Density functional theory
The DFT is a computational method that allows to perform QM modelling of divers systems including condensed-matter physics, quantum chemistry, and the ground-state of a variety of molecules that could be of biological origin. A good review on the subject can be consulted in Ref. [1]. The electron density ρ( r) is the central quantity in DFT, since it determines the probability of finding any of the N electrons of the system within the volume element d r. It is defined as the integral over the spin coordinates of all electrons, together with the normalisation condition, ρ( r)d r = N . The electron density is a positive function that can be experimentally measured and it shows an exponential decay for large distances, ρ( r) ∼ exp[2 √ 2I| r|], where I is the exact ionization energy. The determination of ρ( r) is based on two fundamental theorems established by Hohenberg and Kohn (H-K) in 1964 [2]. For an external potential v( r) acting on the N electrons of the system, the first H-K theorem states that v( r) is determined (to within a constant) by the electron density ρ( r). In other words, v( r) is a unique functional of ρ( r). The second H-K theorem states that for all v-representable densities, being the ground-state energy and ρ 0 is the corresponding ground-state density. The second theorem is of course a direct consequence of the variational principle [3].
Besides the H-K theorems, the Kohn-Sham (K-S) equations can be determined as a function of N orbitals ψ i [4]. Starting with the expression for the energy functional E[ρ] that depends on the interactions between the N electrons and M nuclei of the system, with Z A being the valence of nucleus A, r 12 = | r 2 − r 1 |, and the explicit form of the term E XC [ρ] is unknown in general. Applying the usual constraint ψ i |ψ j = δ ij , leads to the set of equations, where, i are the energy eigenvalues of orbitals ψ i , and V XC = δE XC /δρ. At this point, it is necessary to provide an approximation for E XC . In the generalised gradient approximation (GGA) there is a dependence on densities and their gradients, This term can be combined with the K-S contribution to obtain a hybrid functional, namely, here, a is a fitting parameter. In the K-S scheme, the orbitals are expanded in a set of L predefined basis functions η µ , The Gaussian-type-orbitals are a convenient choice for the basis functions, i.e., η l ∼ r l exp[−κr 2 ]. We obtained the ground-state of a poly-L-lysine (also denoted as polylysine, for short) of 150 residues in vacuum by using the package Orca [5] and the density functional B3LYP [6]. This peptide has two stable configurations, namely, α-helix and β-sheet, which are employed to determine two constants parameters of our CD model. On the other hand, modelling a polylysine in aqueous solution is nowadays out of the scope of DFT not only from the computational point of view. There are still some difficulties in describing certain kinds of intermolecular interactions such as dispersion forces. Instead, we made use of a classical force-field to model our peptides dissolved in water.

Circular dichroism
An important feature of biological molecules is that a considerable number of them are chiral, in other words, they are not identical to their mirror images. Chiral structures can be distinguished and characterised by polarised light [7]. CD is the differential absorption of left and right circularly polarised light or the difference between their respective extinction coefficients, The rotational strength of a transition A ← 0 is the integrated intensity of a single wavelength band in the CD spectrum with µ e and µ m being the electronic and magnetic transition dipole moments. There are theoretical calculations of CD spectra within this QM formalism that are rather involved, but they are only in qualitative agreement with the corresponding experimental data [8] and that is a consequence of the used approximations. In this paper, we present an approach that combines the accuracy of the structures obtained by means of DTF together with experimental data and the results of all-atom (classical) molecular dynamics (MD). The resulting model is computationally tractable even for proteins of hundreds of residues. Fifty years ago, DeVoe proposed a classical theory of CD [9]. For an applied electric field E i , the induced dipole moment µ i at each residue of polarisability tensorα and components α i , is given by where the components of the dipole interaction tensorT are defined as, and (p, q = x, y, z). The formal solution of Eq. (9) isμ =BÊ with µ i = | µ i | e i . It can be assumed that an average polarisability α is able to approximate the three principal components of the monomer polarisabilities [10]. An advantage of the polarisability α is that it can be extracted from experimental data, as explained in the following lines. According to the Beer-Lambert law, the monomer extinction coefficient ε is proportional to the absorbance [7], which is a measurable quantity. Additionally, we use a first order approximation to relate the extinction coefficient with the imaginary part of α, i.e., Here, k 1 is a constant, m is the total number of identical monomers (residues, in our case), and α i = α(ω) with ω being the angular frequency. The real part of α(ω) is obtained through a Kramers-Kronig transform, On the other hand, CD can be determined through the ellipticity [θ] ∼ ∆ε, whereÂ = [α −1 +Ĝ] −1 , G ij = T ij e i · e j , and the constant k 2 has to be determined. The results of this model are presented in the next section.

Experimental part
In this section, we first depict the synthesis of our peptides and the experimental techniques used to describe them. These results together with the DFT and MD calculations are the main ingredients of our model of CD.

Poly-L-lysine synthesis
The monodisperse polymers of L-lysine (of 7, 8, and 20 residues) were manually synthesised by using Fmoc chemistry on 0.22 mmoles of Rink amide 4-Methylbenzhydrylamine hydrochloride copoly (styrene-1% Divinylbenzene) resin (100-200 mesh, loading 0.74 mmol/g, Merck-Millipore). Synthesis was carried out only with L amino acids protected with Nfluorenylmetoxycarbonyl group (Fmoc). The resins were solvated by successive washing A ninhydrin test was used to monitor amino acid couplings and deprotections along the dipeptide synthesis [11]. Finally, the peptides on Rink-MBHA resin were cleaved by treatment with 3 ml of a mixture of Trifluoroacetic acid/water (95:5 v/v) at 0 • C for 15 min and at room temperature for 1.5 hrs. The acid was removed by evaporation and the crude peptide was precipitated with tert-butyl methyl ether/cyclohexane mixture (1:1). The products were characterised by HPLC-MS and ESI mass spectrometry. UV detection was performed at 220 nm, in an Agilent Technologies SQ LC/MSD G6120B instrument. The polymers were purified by HPLC using a Symmetry C18 Preparative Column (100Å, 5 µm, 19 mm × 150, Waters R ) in a Waters instrument with a 2487 Dual λ Absorbance Detector, Waters 600 controller, Waters Delta 600 pump, degasser and Waters Fraction collector III with a manual injector Rheodyne R 3725i controlled by the Empower II software.

Experimental techniques
We obtained a series of UV and CD spectra of our peptides in aqueous solution. The types of spectra were simultaneously acquired using a Jasco J815 spectropolarimeter interfaced with a Peltier unit. Far-UV measurements were performed in 1 mm pathlength quartz cells with PTFE stopper (Hellma). In all cases, three scans were done with a response time of 2 s, a band width of 0.1 nm, and a scan velocity of 20 nm/min.
A range of concentrations between 25-200 µM was used for poly-L-lysine of 7 and 8 residues solution, 10-100 µM for poly-L-lysine of 20 residues and 0.25-2 µM for commercial poly-Llysine. The experiments in aqueous solution were acquired at 295 K, in NaOH (1 mM, pH 10.8) storaged under nitrogen to avoid pH change during the experiments. For each experiment, a blank solution spectrum was subtracted. The spectra were collected between 187-250 nm. The raw ellipticity was converted to molar ellipticity, [θ]. Samples were left 30 min at 278 K to allow equilibration before each experiment.

Results and discussion
The Far-UV spectra are a measure of the absorbance of our peptides in solution, which is equal to the extinction coefficient multiplied by the pathlength and the molar concentration. From the analysed spectra of the three synthesised peptides, we found better results with the polylysine of 7 residues. In the inset of Fig. 1, we present the extinction coefficient ε of this peptide. Assuming an identical polarisability for each residue and using Eq. (11), we obtain the imaginary part of α of a single residue and its real part is computed through Eq. (12). Both components of α are shown in Fig. 1.
In the following stage, the complex polarisability is used together with the structural information of the peptide under study to determine its corresponding ellipticity from Eq. (13). Starting with the optimised polylysines obtained by DFT, we obtain the positions of the centres of polarisability (CP) of each monomer, which are analogous to the centres of mass (instead of using the atomic masses, we employ the dipole polarisabilities of neutral atoms [12]). The electric field E i is null in the z-direction and the vectors e i come from Eq. (9). In order to improve the statistics of [θ], we rotate the peptide along a spiral on the unit sphere, as suggested by Park et al. [13]. The constants k 1 and k 2 are determined by comparing our spectra with the well-established profiles of known structures [7]. As a result, in Fig. 2 we draw the spectra of two polylysines of 150 residues in the configurations of α-helix and β-sheet. As it can be observed in such figure, the two spectra are at least in qualitative agreement with the known spectral shapes. We finally performed a comparison with experimental data. The commercial polylysine has an estimate of 230 residues and we carried out the corresponding MD simulation with the package Gromacs with the Gromos 53A6 force-field and using explicit solvent (three-point SPC/E representation of water molecules). In Fig. 3, we compare the experimental results with the prediction of our model for the molar ellipticity, provided that both systems have similar conditions of pH and temperature. As it can be observed in such figure, there is a fairly good agreement between the modelled spectrum and the experimental one, especially in the region of λ > 210 nm.

Conclusions
We have presented a model of circular dichroism that is able to reproduce, at least in a qualitative way, the known spectra of secondary structures like the α-helix and the β-sheet. In order to calibrate our model, we made use of structures obtained with DFT and experimental data such as the absorptivity and the CD spectra of small peptides that were synthesised by us. These experimental quantities are intended to capture subtle effects that include the interaction between the solvent and the peptide as well as other possible physical interactions (electrostatic, van de Waals, etc.) that still cannot be properly described with QM. More important, the result of our model is in reasonable agreement with an experimental CD spectrum and it was obtained in a reduced computational time.