AMBER-DYES in AMBER: Implementation of fluorophore and linker parameters into AmberTools

Understanding protein dynamics is central to understanding these biomolecules’ functions.^1–3 Investigating protein dynamics on the nanosecond to second or even higher timescales with high resolution is fundamental to detecting conformational changes, folding, association and dissociation of complexes, and enzymatic cycles.^2,4 Molecular spectroscopy is powerful in that sense because it allows observing the required timescales and provides, in addition, a high spatial resolution for molecules in solution.^2,4,5 It is applicable to single-molecule and multi-molecule systems^2,4,6–8 and, therefore, able to investigate a wide variety of biological and biophysical questions governed by protein dynamics. Since naturally occurring fluorescence in the visible range is rare in biomolecules, labels are required. Labels often consist of a fluorescent dye and a flexible linker, covalently connecting the dye to the biomolecule and allowing it to move freely, thereby minimizing quenching mechanisms, e.g., due to surface sticking.^1,9–11

Measuring Förster Resonance Energy Transfer (FRET) is a commonly used fluorescence technique.^12,13 FRET is a nonradiative energy transfer mechanism based on a dipole–dipole interaction between an exited donor dye and a ground-state acceptor dye. The efficiency of the energy transfer (FRET efficiency) is mainly dependent on the distance between the two dyes and on the relative orientation factor κ² of the two electronic transition dipole vectors.^12,13 The FRET efficiency can be determined either by measuring the change in intensity while energy is transferred from the donor to the acceptor or via the change in lifetime of the donor. Precise values of the orientation factor are then needed to convert the FRET efficiency into the distance between the two dyes. That way, FRET can be used as a spectroscopic ruler for quantitative structural studies^1,6,14,15 to measure inter-dye distances in the range of 2 nm–10 nm.^1,2,15 Unfortunately, since detailed information about κ² cannot be measured on the same timescales as FRET, approximations have to be used. Commonly, it is assumed that dyes move isotropically within the donor’s lifetime, resulting in an approximated orientation factor κ²_iso = $23$⁠.^1,2,4 Despite yielding good results for freely moving, non-interacting dyes in solution, this approximation may have shortcomings for protein-bound dyes as such dyes may have a hindered isotropic motion, which may lead to deviations from κ²_iso. Accordingly, the challenge of obtaining detailed information about κ² has been addressed recently by combining experimental data and results from molecular dynamics (MD) simulations.^9,16–18 Explicitly simulating motions of dyes and linkers covalently bound to a biomolecule resulted in accurate anisotropy values, which enhanced the accuracy of inter-dye distance estimates.^10,11,19

A particular challenge for setting up and performing such MD simulations is the availability of an appropriate dye/linker force field that is consistent with the force field used for the protein and compatible with the used simulation software. Graen et al.¹ derived such a force field for commonly used dyes with corresponding linkers that is consistent with the Amber ff99SB force field²⁰ for proteins and its most current successor, the ff14SB force field;²¹ the dye/linker force field has been termed “AMBER-DYES.” Parameters have been made available for 22 commonly used fluorescent dyes of the families Invitrogen Alexa, Lumiprobe Cy, and Atto-Tec Atto, and six linkers (Table I). So far, this parameter set has been available only in the GROMACS^22,23 software suite, which limited its application in the context of the computationally efficient graphics processing unit (GPU) implementation of the AMBER MD engine.^24,25 Here, we report on the adaptation and implementation of the AMBER-DYES parameter set into the free-of-charge AmberTools19 suite,²⁶ which consists of several independently developed packages for biomolecular simulations that work well by themselves and with Amber18.²⁶ The resulting parameter set is termed “AMBER-DYES in AMBER.”

A major challenge for explicitly simulating the dynamics of fluorescent molecules attached to a biomolecule is the system preparation and setup of MD simulations with appropriate force field parameters and topology representations of the molecular components. To address this challenge for MD simulations with AmberTools and Amber, we generated respective parameter and topology files for “AMBER-DYES in AMBER” by adapting the AMBER-DYES parameter set derived by Graen et al.¹ available in GROMACS, as described in the following. Furthermore, to simplify the system preparation and setup of MD simulations with explicit fluorescent molecules, we implemented an easy-to-use workflow.

With the pdb2gmx command of the GROMACS software version 2016.4,^23,27 topology and coordinate files of dyes and linkers (Table I) with AMBER-DYES parameters were created. These files were then transferred into AMBER topology (often called prmtop in AMBER terminology) and coordinate (prmcrd) files using ParmEd version 2.7.2²⁸ from AMBER16.²⁹ Using the AMBER topology and coordinate files, lib files, containing information on connectivity, atom names, atom types, residue names, and charges, and dat files, containing force field parameters for all the bonds, angles, dihedrals, and atom types in the system, were created with ParmEd, which are needed for the general system setup for MD simulations in AMBER applying LEaP²⁶ or other preparatory programs within the AMBER suite. The lib and dat files were modified manually to overcome inconsistencies between the dye and linker representations in GROMACS vs AMBER. As such, the three-letter code for atom names of GROMACS was changed into a two-letter code (Table S1 of the supplementary material) compatible with the AMBER naming scheme. Commonly used force fields such as GAFF,^30,31 GAFF2,³⁰ ff99SB,²⁰ ff14SB,^21,32,33 lipid17,³⁴ and GLYCAM_06j³⁵ can be used together with AMBER-DYES in AMBER. The lib file was modified by defining connecting atoms to create a covalent bond between the dye and the linker. Likewise, the linkers’ base atoms were defined such that they can be covalently connected into the protein chain; note that “linker” here refers to a residue that consists of a natural or non-natural amino acid part to which additional moieties may be connected, e.g., by a reaction of a maleimide with a cysteine (Fig. S1B).

To validate the adapted parameterization, a one-step minimization resulting in essentially a single-point energy calculation of each dye/linker pair in GROMACS or AMBER was performed using ff14SB^21,32 and the AMBER-DYES or AMBER-DYES in AMBER parameters, respectively; for both programs, the same starting structures were used. The minimization was performed in vacuum using an interaction cutoff of 999 Å to include interactions between all atoms of the system. Potential energy terms were then compared after converting GROMACS energies into units used by AMBER. Considering that in GROMACS only integer van der Waals (vdW) scaling factors are allowed, which are used to scale the short-range interactions with respect to experimental data, whereas in AMBER, floating point factors can be used, and that GROMACS and AMBER report energy values with different numbers of decimal places, perfect agreement of total potential energies as well as individual energy components between original AMBER-DYES and adapted AMBER-DYES in AMBER implementations was obtained with a relative error <0.01% for each dye/linker pair (see one example in Table II; for all other pairs, see Tables S2–S127 of the supplementary material).

As pointed out before,³⁶ in the original parameters from the work of Graen et al.,¹ there is an issue regarding the product of the cysteine and the maleimide in every linker ending in a cysteine (namely, C1C, C1R, E1N, C2C, C2R, E2N, C3C, C3R, and E3N linker). Linkers with these original parameters are named with an “^*” at the end to emphasize the presence of the issue. The double bond between the two carbon atoms in the maleimide ring is retained in the product (Fig. S1), although it must become a single bond due to the Michael addition, with two hydrogens bound to the carbon that is not bonded to the sulfur (Fig. S2). These bond properties are neither present in the parameters nor in the PDB structures provided by Graen et al.¹ The parameters and PDB structure provided by Graen et al. lead to a more rigid linker than used in experiments.

To address this issue, we modified the linker structures such that they reflect the correct Michael addition product. The new structures were used to perform a restrained electrostatic potential (RESP³⁷) fitting for atomic charges using Gaussian (Gaussian16³⁸) and Antechamber (AmberTools19³⁹). The capped linkers were refitted, restraining the global charge to zero and all atoms with their respective original charges, except the atoms participating in the obtained single bond and the sulfur (Fig. 1). The newly fitted charges were used to create updated parameters for the linkers (Fig. S2). In the new structures, the bonded sulfur atom does not lie in the ring plane anymore. Thus, the linker is bent closer to the protein surface. Additionally, the carbon bonded to the sulfur is a chirality center. Upon reaction of the maleimide with the cysteine, both configurations occur equally likely; here, we set the configurations of this atom of linkers C1R, C2R, and C3R to R. An overview over representative modified linkers is depicted in Fig. S2. These modified parameters constitute the AMBER-DYES in AMBER implementation. All generated files will be made available within AmberTools19 (http://ambermd.org/AmberTools.php).

LEaP²⁶ is the primary program to create a new simulation system in AMBER. For ease of use, a leaprc.amberdyes file was created, which, when sourced in LEaP, automatically loads the lib and dat files for the respective dye/linker combinations (Table I). With the workflow described in the following, the dye/linker combination is then attached to the protein. Afterward, the generated system, specified in topology (prmtop) and coordinate (prmcrd) files, can be used for MD simulations, as described, for example, at the end.

As an example for the workflow, we use the protein T4 Lysozyme (T4L; PDB ID 148l) together with the Alexa488 (A48) fluorophore and the linker C1R, consisting of a maleimide attached to a cysteine. First, we removed buffer and solvent entries from the PDB file and chose a residue [here, ILE 3; Fig. 2(a)] to be used as the attachment point. To attach the linker, the name of the attachment point residue is changed from “ILE” [Fig. 2(a)] to the linker name “C1R” [Fig. 2(b)], and all atoms of the attachment point residue but the CA atom are removed [Fig. 2(c)]. During system preparation, LEaP will create the missing linker atoms from information in the lib file. The last step of the preparation is to append to the PDB file a new line after the TER card with the residue name of the dye (here, A48) and the atom name “C99” [Fig. 2(d)]. The atom name “C99” for the connecting atom is consistent for all dyes in the AMBER-DYES in AMBER parameter set. Again, LEaP will create the missing dye atoms from information in the lib file.

After the preparation of the PDB file, LEaP²⁶ is used to load the required force fields, solvate the system, and save the files needed for MD simulations. Sourcing leaprc.amberdyes will load all necessary dye/linker parameters [Fig. 2(e)]. Upon loading the solute, in our case, T4 lysozyme with one altered residue “C1R” and the additional “C99” atom of the dye residue A48, missing atoms are automatically added. The connection of the linker and the dye has to be set manually by specifying a bond between the atoms “C99” (dye) and “N99” (linker) [Figs. 2(e) and 2(f)]. The final step is to select the linker and the dye and relax the selected substructure to obtain a feasible dye/linker configuration (Fig. S3A of the supplementary material). Further steps depend on the user’s preferences for setting up MD simulations, such as adding counter ions, selecting the box size, and solvating the solute with a chosen water model.

Following this workflow, the user obtains a protein with an attached dye/linker [Fig. 2(f); Fig. S3A of the supplementary material]. It is possible to attach multiple dyes, also of the same type, by applying the described changes to multiple attachment point residues, appending the PDB file with the corresponding “C99” atoms of the dyes, separated by TER cards, and connecting the linkers and dyes in LEaP. Additional post-translational modifications such as phosphorylation or glycosylation can also be applied due to the orthogonal dye and linker residue and atom names. Although this workflow has worked in all our test cases so far, the user is still advised to carefully check the outcome. Clashes between the protein and the initial placement of the dye/linker are usually resolved by LEaP’s “relax” command or in the first steps of a subsequent energy minimization. In rare cases, however, dye/linker attachments may lead to erroneous configurations (Fig. S3B of the supplementary material), which have to be resolved manually, e.g., by rotating the dye/linker pair to the outside of the system.

As an application example, we again use T4 lysozyme and Alexa488 as a donor dye attached via the C1R linker at residue 19. As an acceptor dye, we use Alexa647 and the corresponding C2R linker attached to residue 119. Following the workflow, we get the structure shown in Fig. 3(a). The structure was solvated in an octahedral box of TIP4P-Ew water with a solvent shell of at least 11 Å around the solute. Na⁺ ions were added at a concentration of 0.15 M. The system was neutralized by adding Cl⁻ ions as counterions and thermalized, and production MD simulations of 2 µs length in the canonical ensemble (constant number of particles, volume, and temperature; NVT) were performed using the particle mesh Ewald method⁴⁰ to treat long-range electrostatic interactions. Hydrogen bonds were constrained using the SHAKE algorithm.⁴¹ An integration time step of 2 fs with a direct-space, nonbonded cutoff of 8 Å was used.

The distance between the centers of the transition dipoles of the two dyes was measured over the course of the trajectory. The locations of the centers of the transition dipoles were defined as in Fig. S4.⁴² The distance distribution was compared to experimental FRET measurements for the same system, which considered the dyes to be freely moving.³ A preliminary assessment of the spatial propensity supports this assumption in this case [Fig. 3(b)]. Accordingly, the computed mean dipole distance of 47.8 Å agrees very well with the experimental mean distance of 49.1 ų obtained from the FRET efficiency and using κ²_iso = $23$ [Fig. 3(c)]. Overall, this shows the applicability of AMBER-DYES in AMBER to investigate protein and dye dynamics within AmberTools and AMBER.

As a further application example, different linker lengths were used to probe their influence on the distribution of dye orientations and potential deviations from κ²_iso = $23$ behavior. The general setup is as before: T4 lysozyme is used as the protein; Alexa488 and Alexa647 are the dyes attached via linkers at residues 19 and 119, respectively; the simulation length is 1 µs each; the TIP4P-Ew^43–45 water model is used.

Comparing the first application example (Alexa488 attached via C1R and Alexa647 attached via C2R) [Fig. 4(a); ⟨κ²⟩= 0.43 ± 0.40; Eq. S1A] to the additional simulations with both dyes attached via the long C3R linker [14 atoms between the backbone carbon to the nitrogen atom bound to the dye; Fig. 4(b)] or the short C2R linker [nine atoms; Fig. 4(c)] reveals differences in the κ² distributions. Using the longer C3R linker [Fig. 4(b)], the distribution of κ² values (computed as described in the supplementary material; Eq. S1A) is shifted to smaller values, and ⟨κ²⟩= 0.35 ± 0.34. This may result from more frequent interactions between the dyes and the protein. In contrast, using the shorter C2R linker [Fig. 4(c)], large κ² values are more frequent, and ⟨κ²⟩= 0.52 ± 0.40, closer to κ²_iso = $23$⁠. Here, the interactions between the dyes and the protein are less frequent.

Note that these results only show the application scope of the AMBER-DYES in AMBER implementation to address questions of the impact of dye motions on the interpretation of FRET experiments. Further simulation efforts with respect to simulation length, numbers of simulations, and different starting configurations are required before general conclusions can be drawn.

As a final application example, different water models were used to investigate their influence on the dye behavior and orientation. The general setup is as before in the first application example: T4 lysozyme is used as the protein, and Alexa488 and Alexa647 are the dyes attached via C1R and C2R linkers at residues 19 and 119, respectively. As to the influence of the water model, with older water models such as TIP3P,⁴⁶ dyes showed a larger tendency to stick compared to experimental data, a problem that is most likely related to the imbalance of protein–water interactions because the dyes share their Lennard-Jones parameters with the protein force field.⁴⁷ Therefore, we now probed the TIP4P-Ew^43–45 and the OPC⁴⁶ water models, which both have lower self-diffusion constants and expansion coefficients,^45,48,49 resulting in a more realistic behavior.

Investigating the time-resolved fluorescence anisotropy (Fig. 5; see the supplementary material for computation, Eq. S1B), only minor differences between the different simulations are visible. The anisotropy decays indicate that both dyes are freely rotating because of the undisturbed anisotropy decay toward zero (Fig. 5) and the good agreement of calculated and experimental distances between the dyes (Fig. 6) using, in the case of experiment, the assumption of κ²_iso = $23$⁠. The anisotropy decays can be fitted as three-exponential decays (Eq. S1C) with correlation times (pre-exponential factors are given in brackets) for Alexa488 of 0.40 ps (0.04), 0.30 ns (0.09), and 1.28 ns (0.27) for TIP4P_Ew and 52 ps (0.04), 0.31 ns (0.07), and 1.28 ns (0.30) for OPC, and for Alexa647 of 33 ps (0.02), 0.37 ns (0.06), and 1.18 ns (0.32) for TIP4P_Ew and 36 ps (0.02), 0.35 ns (0.09), and 1.20 ns (0.29) for OPC. In all four cases, the relaxation times of the MD-derived anisotropy decays are ∼1.2–1.3 ns. These relaxation times are in good agreement with experimental results (1.5 ns) for a free Alexa488 dye.¹⁰ 

Neither of the state-of-the-art water models has a pronounced influence on the distance distribution between the dyes (Fig. 6). Comparing the simulation with TIP4P-Ew water [Fig. 6(a)] to that with OPC water [Fig. 6(b)], both distance distributions are similar. The results of both models are in good agreement to the experimental value of 49.1 Å, indicated by the orange line, and an experimental standard deviation of ∼1 Å, indicated by the orange box.³ Note again that these results only demonstrate the application scope of the AMBER-DYES in AMBER implementation but lack the potential for generalization as to an influence of the water model on dye properties.

In this study, we adapted the AMBER-DYES parameters derived by Graen et. al.¹ into AMBER-DYES in AMBER to generate a force field applicable for commonly used fluorescent dyes and linkers to a protein within the AMBER suite of biomolecular simulation programs, in particular, the computationally efficient GPU implementation of the AMBER MD engine. The adapted parameters result in potential energies, including energy components, that differ by <0.01% from those computed with the original parameters, demonstrating the validity of the adaptation. Following this, we updated the parameters to address for a chemical inaccuracy in the original. Until recently, users of the AMBER suite of biomolecular simulation programs were primarily dependent on self-built solutions or third-party programs to simulate proteins with attached fluorescent dyes. This made the setup of such systems time-consuming and error-prone. As shown in the example cases, the adapted AMBER-DYES in AMBER parameters and the developed LEaP-based workflow now provide a user-friendly and efficient way to generate a wide variety of starting systems, including proteins with multiple identical or different dye/linker combinations. That way, MD simulations of explicit dye representations with AMBER are now easily possible, allowing, e.g., the probing of commonly used approximations for dye localization and/or orientation^2,3 or to model FRET.⁴² 

Although AMBER-DYES in AMBER together with the workflow simplifies a large part of the building process of the fluorescent systems, the user is still required to choose sensible attachment points and critically evaluate if the built system is representative from a structural point of view. Clashes due to dyes pointing inside the protein or the lack of space for the dye around the attachment point, which are not resolved by LEaP’s “relax” command, have to be manually fixed by the user by changing the starting orientation of the dye and/or choosing another sensible attachment point, respectively.

The AMBER-DYES in AMBER force field works together with other commonly used state-of-the-art AMBER force fields such as GAFF, GAFF2, ff99SB, ff14SB, lipid17, and GLYCAM_06j, which allows simulating post-translationally modified proteins and/or protein-ligand complexes and/or proteins in membrane environments. Furthermore, explicit water models such as TIP3P, TIP4P, TIP4P-Ew, and OPC can be used with the force field. Finally, following the parameter development described in Ref. 1 and the adaptation steps described herein, AMBER-DYES in AMBER can be easily extended by further linkers and fluorescent molecules.

See the supplementary material for supplementary methods for computations; supplementary tables of atom names of the dye and linker atoms available within AMBER-DYES and AMBER-DYES in AMBER for AMBER and GROMACS and comparisons of potential energy (components) for AMBER and GROMACS for dye–linker combinations; supplementary texts of input scripts for AMBER and GROMACS for the single-point energy calculations; and supplementary figures showing chemical structures of the linker, worked examples of the workflow, and location of the (centers of the) transition dipoles.

The data that support the findings of this study are openly available in AmberTools19 at https://ambermd.org/AmberTools.php.³⁹ 

We are grateful for computational support and infrastructure provided by the “Zentrum für Informations-und Medientechnologie” (ZIM) at the Heinrich Heine University Düsseldorf and the computing time provided by the John von Neumann Institute for Computing (NIC) to H.G. on the supercomputers JURECA and JUWELS at the Jülich Supercomputing Centre (JSC, user IDs: HDD201 and schepers1). This study was supported in part by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - Projektnummer 267205415 - SFB 1208 (A03 to H.G.).

The authors declare no competing financial interest.