Modeling
The demanding-accuracy-character of the proposed diagnostic project led to the development of models for its full extent, revolving around four main axes: design, analysis, prediction, validation. Firstly, in order to determine the trigger sequence and design the toehold switches, a computational tool was created. Scanning through conserved and unique regions of the MERS-CoV genome and taking into consideration specific thermodynamic criteria, a neural network chose the 4 most promising sequences. Secondly, a kinetics study was performed using appropriate equations, giving insight on the reaction times of the system. Thirdly, the newly proposed reporter system’s (split reporter protein interaction via coiled coils) behaviour had to be predicted via Molecular Dynamics Simulations, and compared to laboratory results. Finally, the development of a final product, a kit that implements all the processes, had to be meticulously studied. Through PDEs, optimum times were calculated ensuring optimum conditions for every step.
Toehold Design
Toehold Parts
-
TBS: Trigger Binding Site, a 30 nucleotides long sequence, reverse complementary to the target RNA. Target RNA will be attracted and finally bind to it.
-
RBS: Ribosome Binding Site, a 12 nucleotides long sequence, responsible for the recruitment of a ribosome during the initiation of protein translation.
-
Linker: A 21 nucleotides long sequence, connecting the reporter protein to rest of the toehold, ensuring minimum interaction between the two.
Target RNA search
The whole Toehold Switch mechanism aims at the recognition of a specific RNA sequence, the target RNA. This 30nt long sequence is part of the MERS-CoV genome. Every target RNA corresponds to a single Toehold. All the possible pairs are created and tested based on thermodynamic criteria. The first target RNA corresponds to the first 30 nucleotides, the second to nucleotides 2-31 etc. That means that practically for a genome of length N, there are N possible switches (pair of target and corresponding toehold). MERS-CoV has a length of about 30.000 nucleotides, rendering this not a negligible computational problem.
Conserved Regions
Different strains of viruses may present significant differences on various regions of their RNA sequence. Nevertheless, conserved regions exist and in fact are the ones targeted in classic diagnostic molecular methods, such as PCR. Consequently, search for target RNA is limited to conserved sequences and not the whole MERS-CoV genome reducing the number of switches tested to about 1000 from 30.000. The computational load is significantly decreased but so does the possibility of finding an admissible sequence. The conserved regions are determined through bibliography[1].
Toehold specificity to MERS-CoV
Due to the relatively large size of the target RNA, the possibility of the exact same sequence to exist in other viruses or humans is relative small. Therefore, BLAST search is conducted only at the set of the most promising targets for time efficiency, confirming the presence of the sequence only at the MERS-CoV genome.
Thermodynamic Evaluation
Every switch needs thermodynamic evaluation to determine its ability to stay closed in absence of target and to open when the target is present. After the evaluation, every switch is scored and the 4 most promising where chosen for in vitro testing.
MFE
MFE(Minimum Free Energy) criterion, ensures the unavailability of the RBS in absence of the target RNA. Highly negative values of MFE imply great stability of the toehold thus zero expression of the reporter protein. On the contrary, values near zero indicate expression leakage, an undesirable phenomenon where reporter protein is expressed in absence of the target, resulting in false positive diagnosis.
ΔMFE
ΔMFE criterion, ensures the binding of the target to the TBS. Two possible states exist for the system containing the target and the toehold, them binding together to form a duplex or them existing separately. ΔMFE will determine the most favorable state, in terms of thermodynamic stability. \( \Delta MFE = MFE_{target} + MFE_{toehold} – MFE_{duplex}\), the higher its value the more possible it is for the target to bind to the TBS.
ΔG
ΔG criterion, ensures the binding of the ribosome to the RBS. Even after the initial unfolding of the switch with the binding of the target to the TBS, the hairpin at the top of the switch may remain closed, detering the binding of the ribosome. High values of ΔG (close to zero, as ΔG is a negative quantity) suggest easy unfolding of the upper hairpin and thus, RBS availability.
Switch-Switch/Target-Target
During actual implementation, the sample will contain a large number of switches. Pairs of triggers-triggers and switches-switches will be created, decreasing, possibly in a significant manner, the expression of the reporter protein. The more negative the MFE of these pairs is, the more possible this unfavorable state is.
Scoring Method
After implementation of the thermodynamic criteria, every toehold corresponds to a set of 5 energies. Those 5-dimensional data can’t give a definitive answer as to which switch is preferable, because no switch will have the most desirable values for every criterion. A score method needs to be introduced, in order to quantify the suitability of every switch. A neural network was trained to accomplish this task.
Neural Network Training
The training of a neural network requires a big amount of data. Experimental data on toehold switches’ performance are available online from Green[2]. Using these data, a neural network was trained with the Levenberg-Marquardt backpropagation method to correlate a set of 5 energies to a score. One layer with 5 neurons is used.
Neural Network Efficacy
The performance of the NN on the training set is more than sufficient. Nevertheless, tests must be conducted to ensure that the performance is satisfying on other data sets too. For this reason 30% of the available data weren’t used for training but as independent test, to ensure that the NN is not overfitted.
Kinetic Modeling
In order to gain a better understanding of the RNA-induced glucose production occurring in the cell-free system in use, we developed a simple kinetic model. In that way, we can use our insight on the detection mechanism to assist our suggested kit design. We also tried to incorporate the toehold in-silico thermodynamic evaluation to the kinetic performance of the riboregulator. We formulated a set of ODEs describing our biochemical reactions and solved them with MATLAB.
Our Aim
By applying the kinetic model that we developed, we wanted to better understand the in vitro detection mechanism itself and the course of our wet lab experiments, as well as assist the kit design process.
Brief Model Description
In our system, the trigger RNA sequence concentration [Trigger] works as the input, while the glucose levels [Glucose] as the output. The model accounts for the DNA circuit transcription (TX) to toehold-regulated trehalose-encoding mRNA, the mRNA binding to the trigger sequence, the trehalase translation (TL), the trehalose enzymatic hydrolysis towards glucose and the RNA denaturation. As suggested in the existing literature on PURE TX-TL system modeling, we assign Michaelis-Menten kinetics to the TX and TL reactions occurring, in order to take into account the observed saturation in the reaction rate[3].
Since trehalase is reported to have an optimum pH of 5,5 and shows lower to no activity in neutral pH, we decided to split our study in two parts; during the first part the trigger is added to the cell-free system where toehold transcription, toehold-trigger binding and trehalase translation take place, while in the second part the produced trehalase is activated by lowering the pH and the glucose production occurs. Schematically, the reactions taking place are:
1st Part
2nd Part
Transcription reactions (1)
The circuit DNA TX follows Michaelis-Menten kinetics. The constants’ values used are found in literature[3, 4]. The suggested \( k_1^{TL}\) value for GFP encoding DNA transcription is 18,2nM/min[4]. Since the trehalase encoding gene is approximately 2,2 –fold bigger in size than the one encoding GFP, we assumed that the TX speed is accordingly slower. The \( K_m^{TX}\) value used [3], is calculated in a PURE system for GFP transcription. Using the same system for our experiments and assuming that there is no strict correlation between the gene type and the transcription enzymes affinity, we can use the same value for TreA-gene transcription.
Translation reactions (3) and (4)
Similarly, the RNA translation step follows Michaelis-Menten Kinetics. We assume that the activated dimer’s translation constant \( k_2^{TL}\) is 2.2-fold smaller than the corresponding constant for the GFP mRNA2. The toehold TL constant is initially estimated to be \( k_1^{TL} \approx \frac{k_2^{TL}}{40}\). We take the assumption that the TL Michaelis-Menten constant corresponds to the ribosome binding affinity to the mRNA molecule, thus is independent to the protein type. The constant also has the same value for the activated dimer and the toehold[3]: \( k_{m_1}^{TL} = k_{m_2}^{TL} = k_m^{TL} \)
Activation reaction (2)
We suppose that the toehold- trigger activation reaction quickly reaches the equilibrium state[3]:
Where \( K_d \) is the activated dimer’s dissociation constant. Its value is thermodynamically estimated a priori:
Where the \(\Delta G^o\) value, corresponds to the difference in the Minimum Free Energy values of the molecules in the two states and in 25 oC[5]:
Initially, we believed that such an assumption could provide us a reliable value for the dissociation constant. Dr Estevez-Torres, suggested that we revise that assumption, because despite its simplicity it might be proven insufficient. For that reason, we compared values thermodynamically derived for toehold/trigger pairs that have been tested in PURE TX-TL systems by Senoussi et.al. with the experimentally derived values[3]. Since the differences were substantial, we decided to initiate a correction multiplication factor that links the thermodynamic estimates to reasonable kinetic constants:
The RNA molecules simultaneously co-exist in many different possible conformations. This explains why the Minimum Free Energy value cannot be used alone for the equilibrium constant estimation.
Denaturation/Deactivation reactions (5) and (6)
We assume that only the RNA species are susceptible to denaturation during the detection test duration. The mRNA follows a first order denaturation kinetics[4].
Trehalose Hydrolysis (7)
Trehalose enzymatic hydrolysis to glucose follows Michaelis-Menten kinetics. In order to perform our calculations, we used the Michaelis-Menten kinetic constants, as calculated by Gibson et.al [6]:
| Constant | Value |
|---|---|
| \( k_1^{TX}\) | 0.138 nM/s |
| \( K_m^{TX}\) | 4.2 nM |
| \( k_1^{TL}\) | 0.003 nM/s |
| \( k_2^{TL}\) | 0.122 nM/s |
| \( k_m^{TL}\) | 265 nM |
| \( R \) | 1.985∙10-3 kcal∙mol-1∙K-1 |
| \( T \) | 310 K |
| \( a \) | 0.11 |
| \( k_d \) | 7.8∙10-4 min-1 |
| \( K_{cat} \) | 199-1 |
| \( K_m \) | 0.41 nM |
The in vitro mechanism
For the in vitro toehold testing we assume that the trigger RNA and the toehold DNA sequences are simultaneously added to the system (time = 0 min). In the following diagrams we present the cell free translation-transcription reactions’ propagation for different Trigger concentrations. After a time course of 60 minutes, we assume an instant pH drop from 7 to 5,5 whereas the trehalose hydrolysis occurs .
Kit Design
The kinetic modeling of the kit set as a primary goal the definition of the time spans of each one of the test steps. Since we did not have enough data to model the kinetics of the Recombinase Polymerase Amplification (RPA), our study focused on the post-amplification steps. We assumed that after the amplification the expected Trigger sequence concentration may vary between 10 nM - 10 μΜ. The amplification step is set to 30 minutes, as to ensure that the input trigger concentration for the cell-free reactions will lie between the above-mentioned range.
After a series of meetings (Integrated HP) we concluded on setting the 1 hour as the upper limit for the whole detection process, since lengthier applications may be proven unsuitable for PoC applications.
As Professor Keith Pardee suggested, we run the amplification and the DNA-circuit-to-toehold transcription in different compartments because they are shown to be incompatible. In order to save time, we run these reactions in parallel, supposing that they instantly mix after 30 minutes. The initial Toehold-DNA concentration in all the cases is assumed to be 3 nM.
Both the TX-TL reaction and the Enzymatic Hydrolysis have to take place in a total of 30 minutes. In order to decide the division between the two steps, we expressed the final glucose levels as a function of the TX-TL step, which produces the hydrolytic enzyme- Trehalase.
During this approach, we set the value of 4 mmol/L glucose as the viral detection threshold, since this value corresponds to typical fasting levels in blood for a healthy adult, targeting to such a value, we ensure that it will be displayed by the glucometer regardless its type and detection range.
In order to reach final glucose levels above the aforementioned glucose limit, the duration of the TX-TL step should lie between 8,8 min - 21 min. While our first pick was the 15 min duration, where the glucose levels reach a pick, we chose to set the TX-TL reaction in 20 minutes, in order to take into account possible negative feedback because of the mixing with the amplification reactants. It is worth mentioning, that the use of round time fractions for each step is essential in order to facilitate the kit operation.
With that being said, The TX-TL step is set in 20 minutes and the enzymatic reaction step in 10 minutes:
3D Simulations
Protein Complex Simulation
Our main goal for performing the Molecular Dynamics simulation of the Split Trehalase Reporter System was to develop a model in order to observe the protein interaction in a solution and, thus, have an insight on the system’s behaviour.
In order to simulate how the two parts of trehalase assemble, when each of them is merged with half a leucine zipper, we first needed to create a model for both parts. Trehalase’s three-dimensional structure, has been crystallographically determined and the leucine zipper is quite a common structural motif in proteins, so it was easy to create models using HHpred and MODELLER[7].
HHpred is a fast server for remote protein homology detection and structure prediction and implements pairwise comparison of profile hidden Markov models (HMMs) [8]. HHpred’s results can be easily forwarded to MODELLER, which is, as the name would suggest, used for homology or comparative modeling of protein three-dimensional structures [9].
Below, you can find listed the sequences that were used for the creation of the models, the PDB structures that HHpred proposed and we selected for the creation of the models and the respective info and statistics:
-
TreA-a part (translated by the toehold switch), contains: ATG + Toehold Linker + Zipper + Zipper Linker + TreA-a (see Design section for more information.
-
HHpred selected results:
Figure 12: On the top, the sequence of the junk amino acids and the leucine zipper part are depicted in grey, while the trehalase part is depicted in blue. Below, the PDB structures are depicted in red or orange. The length of each line indicates the part of the sequence that each structure covers on the model. In the middle of each line, the structure’s PDB ID is indicated with white letters.
Figure 13: Predicted structure of TreA-a part along with its zipper, the long purple a-helical part. Visualization enabled via the VMD package. -
-
TreA-b part contains: Zipper + Zipper Linker + TreA-b (see Design section for more information)
-
HHpred selected results:
Figure 14: On the top, the input sequence is depicted in blue. Below, the PDB structures are depicted in red. The length of each line indicates the part of the sequence that each structure covers on the model. In the middle of each line, the structure’s PDB ID is indicated with white letters.
-
A molecular dynamics (MD) simulation for the aforementioned complex of the Split Trehalase Reporter System was performed using the GROMACS packageup [10][11], taking into account the MD simulations performed by Aronsson, et al. 2012 [12]. The stages of the simulation, observations as well as comments are mentioned in the “Information & Comments” section below.
Note 1: The NAMD package could also be used, if one is more familiar with it. However, GROMACS is preferred for protein simulations, as it displays better performance.
Note 2: Because of the size of the system and the amount of time needed to produce promising results, the whole process run on the national supercomputer, ARIS of the Greek Research & Technology Network (GRNET).
Information & Comments:
-
In order to facilitate the process, the two protein parts were brought to an 8-12Å distance from each other, with appropriate orientation. This was mainly done so that the size of the water box (discussed below) would be manageable. This means that the solvated system would contain a reasonable number of atoms that the computer could process in a relatively rational amount of time. At the same time, it should be mentioned that the position and orientation of the proteins could, but it is not necessary, aid their interaction.
-
Force Field: The CHARMM36 [13][14] force field was implemented for the peptide and ions parameters.
-
Water box, Solvation, Charge: A rhombic dodecahedron box was chosen, with a periodic distance equal to the one of a cubic box but with ~71% its volume, ensuring minimal number of atoms of the solvated system. The three-point TIP3P [15] water model was used. The total charge of the system was set to neutral by adding Natrium ions.
-
Other Details: Periodic boundary conditions (pbc = xyz) were implemented on all stages. Regarding the electrostatics, the particle-mesh-Ewald method was chosen, setting the grid length ( fourierspacing) equal to 0.16nm, with a 1.2nm cut-off. A cut-off of 1nm was included for the Van der Waals interactions.
-
Minimization: A steepest descent algorithm was used for energy minimization, with a maximum force tolerance (emtol) of 1000.0kJ mol-1 nm-1. As shown on the diagram below, the system’s energy is minimized in less than 1.5ns.
Figure 16: Potential energy vs. Time for the duration of the Energy Minimization process. Plot obtained via XMGrace. -
Temperature Equilibration (NVT): The temperature of the peptides and solvent was kept constant at 310K, at the temperature our real-life system is expected to work, with V-rescale as the chosen thermostat. Each group of molecules (Protein, Water_and_ions) was coupled independently. The temperature time constant was 0.1ps. Position restraints were implemented for the protein, while the LINCS constraints algorithm was used for h-bonds constraints.The integration of the equations of motion was performed using a lead frog algorithm (md) with a time step of 2fs, for a total time of 10ns.
Figure 17: Temperature vs. Time for the duration of the NVT Equilibration process. Plot obtained via XMGrace. -
Pressure Equilibration (NPT): Temperature coupling was on during this stage, using V-rescale thermostat, as described above. For pressure coupling, the Berendsen barostat was chosen, which is best for the equilibration stages. The reference pressure was 1.0bar, the pressure time constant was 2ps, while the compressibility was 4.5e-5 (appropriate value for T=310K or T=300K and p=1.0bar). Position restraints were implemented for the protein, while the LINCS constraints algorithm was used for h-bonds constraints. The integration of the equations of motion was performed using a lead frog algorithm (md) with a time step of 2fs, for a total time of 50ns. Suggestion: A suggestion our team received was to pressure equilibrate in two steps, one with position restraints and one without, and even reducing the total NPT simulation time down to 10ns. The suggestion was tested for reference on a different system, showing good results.
Figure 18:Pressure vs. Time for the duration of the NVT Equilibration process. The fluctuations of pressure are to be expected, with the running average fluctuates to a less degree. Plot obtained via XMGrace. 
Figure 19: Density vs. Time for the duration of the NVT Equilibration process. The density converges nicely, with little fluctuations that are to be expected. Plot obtained via XMGrace. -
Production Run Temperature and Pressure coupling were on during this stage, with a V-rescale thermostat and a Parrinello-Rahman barostat. The system was set free of any position restraints. The LINCS constraints algorithm was used for h-bonds constraints. The integration of the equations of motion was performed using a lead frog algorithm (md) with a time step of 2fs, for a total time of 200ns.
-
Analysis: Below the trajectory of the protein complex is showed, for a total duration of 100ns (the last 100ns of the simulation).
One of the first thing one notices in the movie above is that the protein molecules drift away from each other and exit the water box. However, this is not entirely true. As mentioned, Periodic Boundary Conditions (pbc) were implemented in the system. That means that when a protein “exits” from the right, its periodic image will “enter” from the left, as protein molecules are free to diffuse in the assigned waterbox. For example:
So, by enabling periodic boundary conditions on the visualization software, one can observe that the proteins actually do approach each other. However, what we have not observed yet, is the leucine zippers interaction.
Observations
As expected, the proteins move around in the water box. Thus, their interaction relies mostly on diffusion terms, on their structure, their polarity and size. In order to observe a significant result, the Production Run stage of the simulation has to run for a significant amount of time, where a 500ns total duration is estimated to produce better results. However, for the total time of the simulation, one can observe that the two molecules do approach each other, then drift apart, then approach again with different orientations, indicating the expected interaction.
Future Plans
-
In order to better observe the interaction of the two leucine zippers and the protein parts, a Steered Molecular Dynamics Simulation could be performed. Through the Steered Molecular Dynamics Simulation, we could also obtain a better estimate of the Gibbs Free Energy and subsequently the equilibrium constant Keq.
-
Using the Protein-Protein Docking approach prior to the MD Simulation, via FRODOCK [16] or HADDOCK [17] servers, a comparison to the method described above can be made. This way, one can determine whether docking is a promising preparatory method for a protein complex prior to the MD Simulation.
-
Due to software limitations, RNA 3D Simulations are very difficult to perform. It is our team’s desire to develop a 3D model for our RNA toehold switches, when the necessary softwares become available.
Temperature Simulations
For the steps of viral RNA amplification and transcription-translation of the cell-free system, a constant temperature of 37oC needs to be achieved. Positive Temperature Coefficient (PTC) thermistors are going to be used for that purpose. Simulations need to be conducted in order to determine the time needed for the kit temperature to reach steady state conditions. A 3D model was formulated, and time dependent simulations revealed the requested time
Temperature is monitored in the center of the large paper disc. The results shown on the image below are used to calculate the steady-state time. It is obvious that the temperature profile is steady after 10-15 minutes, thus 10 minutes are deemed enough as a waiting time before the start of step 3 (RNA amplification). Furthermore, using the same diagram we investigate the solution dependency on the used mesh. The blue line (Original Mesh, 18.000 elements) and the red circles (Finer Mesh, 45.000 elements) completely overlap, so we deem that the normal mesh is sufficient enough and more practical, thus it is going to be used for the further calculations needed.
Heater position is of great importance for the development of the temperature profile. A parametric study was conducted to obtain the most favorable, in terms of fast achievement of steady-state conditions, position relevant to the initially suggested. All the results presented used this optimum position.
Humidity of the paper may play a significant role at the heat transfer inside the paper and consequently at the time needed for reaching steady-state conditions. Humidity affects the heat transfer coefficient (k), but the correlation of the two can’t be described adequately due to the bizarre nature of the liquid. Instead we conducted a parametric investigation with respect to the heat transfer coefficient directly. The results obtained clearly indicate that the 10 minutes calculated in the beginning are enough no matter the value of the heat transfer coefficient.
Based on the above investigation, a time of 10 minutes is needed after step 2 in order to achieve a constant temperature of 37oC, required in step 3 (RNA amplification). Step 4 (transcription-translation of the cell-free system) also requires a constant temperature of 37oC, but this time the temperature is already achieved so zero lag time is required.
References
[2]Green AA, Silver PA, Collins JJ & Yin P., “Toehold switches: De-novo-designed regulators of gene expression”, Cell, 159(4):925–39, (2014).
[3]Senoussi, A. et al. Quantitative characterization of translational riboregulators using an in vitro transcription-translation system. bioRxiv 290403 (2018). doi:10.1101/290403.
[4]Stögbauer, T., Windhager, L., Zimmer, R. & Rädler, J. O. Experiment and mathematical modeling of gene expression dynamics in a cell-free system. Integr. Biol. 4, 494 (2012).
[5]Zadeh, J. N. et al. NUPACK: Analysis and design of nucleic acid systems. J. Comput. Chem. 32, 170–173 (2011).
[6]Gibson, R. P. et al. Molecular Basis for Trehalase Inhibition Revealed by the Structure of Trehalase in Complex with Potent Inhibitors. Angew. Chemie 119, 4193–4197 (2007).
[7] Zimmermann, L., Stephens, A., Nam, S.-Z., Rau, D., Kübler, J., Lozajic, M., … Alva, V. (2018). A Completely Reimplemented MPI Bioinformatics Toolkit with a New HHpred Server at its Core. Journal of Molecular Biology, 430(15), 2237–2243. doi:10.1016/j.jmb.2017.12.007.
[8] Soding, J., Biegert, A., & Lupas, A. N. (2005). The HHpred interactive server for protein homology detection and structure prediction. Nucleic Acids Research, 33(Web Server), W244–W248. doi:10.1093/nar/gki408.
[9] Webb, B. and Sali, A. 2016. Comparative protein structure modeling using MODELLER. Curr. Protoc. Protein Sci. 86:2.9.1-2.9.37. doi: 10.1002/cpps.20
[10] Hess, B., Kutzner, C., van der Spoel, D., Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 4(3):435–447, 2008.
[11] M.J. Abraham, D. van der Spoel, E. Lindahl, B. Hess, and the GROMACS development team, GROMACS User Manual version 2018.3, www.gromacs.org (2018)
[12] Aronsson, C., et.al. Self-sorting heterodimeric coiled coil peptides with defined and tuneable self-assembly properties. Scientific Reports, 5. doi:10.1038/srep14063 (2015).
[13] Huang, J., Rauscher, S., Nawrocki, G., Ran, T., Feig, M., de Groot, B. L., … MacKerell, A. D.. CHARMM36m: an improved force field for folded and intrinsically disordered proteins. Nature Methods, 14(1), 71–73. doi:10.1038/nmeth.4067 (2016).
[14] Best, R. B., Zhu, X., Shim, J., Lopes, P. E. M., Mittal, J., Feig, M., & MacKerell, A. D.. Optimization of the Additive CHARMM All-Atom Protein Force Field Targeting Improved Sampling of the Backbone ϕ, ψ and Side-Chain χ1 and χ2 Dihedral Angles. Journal of Chemical Theory and Computation, 8(9), 3257–3273. doi:10.1021/ct300400x (2012).
[15] Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W., & Klein, M. L. . Comparison of simple potential functions for simulating liquid water. The Journal of Chemical Physics, 79(2), 926–935. doi:10.1063/1.445869 (1983).
[16] E. Ramírez-Aportela, J.R. López-Blanco, and P. Chacón (2016). FRODOCK 2.0: Fast Protein-Protein docking server. Bioinformatics, 32(15), 2386-2388.
[17] Cyril Dominguez, Rolf Boelens and Alexandre M.J.J. Bonvin. HADDOCK: a protein-protein docking approach based on biochemical and/or biophysical information.J. Am. Chem. Soc. 125, 1731-1737 (2003).
