Front. Mol. Biosci. Frontiers in Molecular Biosciences Front. Mol. Biosci. 2296-889X Frontiers Media S.A. 10.3389/fmolb.2021.659687 Molecular Biosciences Original Research Molecular Mechanics Study of Flow and Surface Influence in Ligand–Protein Association Kaushik Shivansh Chang Chia-en A. * Department of Chemistry, University of Chemistry, Riverside, CA, United States

Edited by: Joanna Trylska, University of Warsaw, Poland

Reviewed by: Nanjie Deng, Pace University, United States; Peter M. Kekenes-Huskey, University of Kentucky, United States; Sanjib Senapati, Indian Institute of Technology Madras, India

 *Correspondence: Chia-en A. Chang chiaenc@ucr.edu

This article was submitted to Biological Modeling and Simulation, a section of the journal Frontiers in Molecular Biosciences

 10 05 2021 2021 8 659687 28 01 2021 06 04 2021 Copyright © 2021 Kaushik and Chang. 2021 Kaushik and Chang

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Ligand–protein association is the first and critical step for many biological and chemical processes. This study investigated the molecular association processes under different environments. In biology, cells have different compartments where ligand–protein binding may occur on a membrane. In experiments involving ligand–protein binding, such as the surface plasmon resonance and continuous flow biosynthesis, a substrate flow and surface are required in experimental settings. As compared with a simple binding condition, which includes only the ligand, protein, and solvent, the association rate and processes may be affected by additional ligand transporting forces and other intermolecular interactions between the ligand and environmental objects. We evaluated these environmental factors by using a ligand xk263 binding to HIV protease (HIVp) with atomistic details. Using Brownian dynamics simulations, we modeled xk263 and HIVp association time and probability when a system has xk263 diffusion flux and a non-polar self-assembled monolayer surface. We also examined different protein orientations and accessible surfaces for xk263. To allow xk263 to access to the dimer interface of immobilized HIVp, we simulated the system by placing the protein 20Å above the surface because immobilizing HIVp on a surface prevented xk263 from contacting with the interface. The non-specific interactions increased the binding probability while the association time remained unchanged. When the xk263 diffusion flux increased, the effective xk263 concentration around HIVp, xk263–HIVp association time and binding probability decreased non-linearly regardless of interacting with the self-assembled monolayer surface or not. The work sheds light on the effects of the solvent flow and surface environment on ligand–protein associations and provides a perspective on experimental design.

 molecular modeling molecular recognition drug design GeomBD ligand-receptor binding GM-109045 MCB-1932984 National Institutes of Health10.13039/100000002 National Science Foundation10.13039/100000001 Introduction

Molecular association is the first critical step in all chemical and biological processes such as the immune response, signal transduction, drug-protein binding, and chemical catalysis (Ozbabacan et al., 2010; Dill and Bromberg, 2012; Baron and McCammon, 2013; Lin et al., 2020). Simulation techniques play a crucial role in investigating the environment that may affect ligand–protein binding, which provides a fundamental understanding of molecular recognition and reduces the cost and time in molecular design. Although diffusion-controlled association rate constants may be approximated analytically, most ligand–protein systems have slower association rates than the diffusion-limited rate because the association event involves multiple steps (Di Cera, 2017; Pang and Zhou, 2017). Conformational rearrangement of both molecules largely determines their binding kinetics, but the two molecules must have an initial encounter first. This first step may be greatly affected by the environment, which can result in different measured association-rate constants.

Although many ligand–protein bindings take place in a closed system without water flow such as in an experimental beaker or cell, ligand–protein association can also occur in a more dynamic environment. Different compartments within a cell also create various membrane environments when ligands and proteins associate and function (Zotter et al., 2017). For example, techniques using continuous flow biocatalysis have been developed recently to improve the efficiency of chemical synthesis, such as improved mixing, mass transfer, and automation (Planchestainer et al., 2017; Britton et al., 2018). Surface plasmon resonance (SPR), a powerful technique to measure molecular binding kinetics and thermodynamics, also utilizes flow chemistry and a continuous flow environment (Hinman et al., 2018; Prabowo et al., 2018; Ershov et al., 2020). However, how the flow may affect the ligand–protein association processes is unclear. Moreover, these methods need to immobilize one molecule on a surface that may have intermolecular interactions with the molecules to be tested. Therefore, the choice of surface and flow rate usually need to be optimized for various systems.

In addition to experimentally measured values such as association rate constants, molecular simulation can reveal atomistic details of molecular association to further interpret experimental results and understand binding mechanisms. Molecular encounter usually involves two molecules searching for each other in a vast space and for longer than a nanosecond search time. Brownian dynamics (BD) simulations have been used for computational assessment of the encounter processes for several decades (Northrup et al., 1984; Huber and McCammon, 2019). Many software packages are available for studying a variety of problems related to bimolecular association. For example, MacroDox, UHBD, BrownDye Simulation of Diffusional Association (SDA), BD_BOX, BROMOC, ReaDDy, Smoldyn, SEEKR, and BDpack are used to probe ligand–protein associations with a flexible or rigid biomolecular model (Madura et al., 1995; Northrup et al., 1999; Huber and McCammon, 2010; Długosz et al., 2011; Schöneberg and Noé, 2013; De Biase et al., 2015; Martinez et al., 2015; Saadat and Khomami, 2015; Andrews, 2017; Votapka et al., 2017). Our group has been developing the GeomBD program, which focuses on using BD to investigate inter-enzyme intermediate transfer and substrate association on surface environments in biosensor and nanoenzyme structures (Roberts and Chang, 2015, 2016). Because of the diverse bio-systems and binding environments, modifying an existing BD package to answer various questions is common practice.

Here we used the HIV protease (HIVp) and xk263 as a model system to study the effect of flow and surface on ligand–protein associations. HIVp belongs to the class of aspartyl proteases and is essential for maturation and assembly of infectious virions (Kohl et al., 1988). The protein cleaves the large polyprotein precursors to mature viral proteins (Huang et al., 2014), an essential function for viral replication, and is a major drug target for AIDS treatment. This is a well-studied system with rich experimental binding data for many inhibitors and FDA approved drugs (Ghosh et al., 2016). The flexible flaps of HIVp can serve as a gate: its opening and closing affects ligand binding. In proteins, an open/closed rate of a gate related to the diffusion of their binding partners determines a fast or slow gating for ligand binding (Szabo et al., 1982; McCammon, 2011). Because of the complex gating behavior of HIVp and its importance in therapeutics, earlier work applied BD to study drug–HIVp binding and also developed a specialized coarse-grained model for HIVp to model the large-scale motions of the flaps (Tozzini and McCammon, 2005; Chang et al., 2006; Kang et al., 2011; Li et al., 2012; Bernetti et al., 2019).

In this work, we used BD simulations to examine the xk263–HIVp association under the influence of ligand diffusion flux and a surface environment. The simulation applied rigid-body BD movements and considered atomistic details using pre-computed grids when computing intermolecular interactions and driving forces. The HIVp was immobilized on a surface or artificially placed 20 Å above the surface. We introduced a self-assembled monolayer (SAM) with a CH3 terminal group to model a hydrophobic surface environment (Cholko et al., 2019). We analyzed ligand association time and binding probability with different diffusion fluxes of xk263, surface and HIVp orientations. The non-polar SAM provided only weak intermolecular attractions with xk263, but the surface did not accelerate the xk263 encounter processes. The x-direction diffusion flux of xk263 significantly reduced weak intermolecular interactions and searching near the surface which shortened the association time but also significantly reduced the probability of successful binding. The concentration gradient of xk263 was affected by the xk263 diffusion flux as well. Our studies suggest that the flow may have noticeable effects on molecular encounter processes and provide insights into the differences in the ligand–protein association in diverse conditions, such as in static cells or a continuous flow biocatalysis environment.

 Methods Model System

The model system is a rectangular prism (400 × 400 × 220 Å3), closed at the top, and consisting of HIVp, xk263, and a SAM surface (Figure 1A). The crystal structure of HIVp and the xk263 inhibitor were obtained from the Protein Data Bank (codes 1HHP and 1HVR, respectively) (Spinelli et al., 1991; Lam et al., 1994). The HIVp flaps, two polypeptides that cover the active site, are in the open position. The semi-open form crystal structure of HIVp was simulated to derive the open form (Supplementary Methods) (Huang et al., 2017). HIVp was kept fixed in the model (Figure 1B). The SAM surface with HIVp at the center consisted of undecanethiol chains on a gold sheet of 1 atom thickness (Supplementary Methods) (Cholko et al., 2019). The size of the SAM was 400 × 400 Å2 having a hexagonal packing pattern with a packing density in the order of 1014 cm−2, which gives an average chain separation of 4.98 Å. The system had a periodicity in the y-direction and termination boundary in the +x-direction as the ligand flows in the +x-direction. HIVp had two orientations (Figure 2).

Model system. (A) HIVp is placed at the center of the CH3-SAM surface with xk263 ligands in a yz-plane when a simulation starts. (B) Structure of an open-form HIVp in cartoon representation. (C) xk263 ligand, and (D) undecanethiol molecule used to build the SAM surface. Notably, in each BD run, the molecular system has only one xk263 and one HIVp. To speed up the calculations, multiple replicas of xk263 are simulated simultaneously, but the replicas do not interact with each other.

Two orientations of HIVp. The protein is placed 20 Å above the SAM (A) parallel to the yz-plane and (B) perpendicular to the yz-plane.

 Simulation Details

The in-house modified GeomBD2 program was used for all BD runs (Roberts and Chang, 2016). The program first creates three grid files: the exclusion volume, screened Coulumbic potential, and 12-6 Lennard-Jones potential for HIV and SAM. Intermolecular interactions between ligand xk263 and HIVp or SAM are grid-based calculations, and the forces are estimated numerically (Roberts and Chang, 2016). To speed up the calculations, these precalculated grids are extended up to 40 and 15 Å around HIVp or SAM for screened Coulumbic and 12-6 LJ potentials, respectively. The force field parameters were obtained from Amber-ff14SB for HIVp and CH3-SAM. The partial charges of xk263 and SAM were computed by using the AM1-BCC charge method with the antechamber program (Cholko et al., 2019). The motion of the ligand is a Brownian motion, governed by an overdamped Langevin equation, in implicit water (Northrup et al., 1984).

ri(t+Δt)=ri(t)-DikBT∂E∂riΔt+2DiΔtR

where Di is the translational or rotational diffusion coefficient, kB is Boltzmann's constant, T is temperature, δE/δri is the potential gradient computed numerically based on the potential grids, Δt is the time step, and R is the stationary Gaussian random number with a zero mean. To obtain the flow for the ligand, the overdamped Langevin equation is modified by adding N, a small fractional number, in the displacement of the x-direction.

xi(t+Δt)=xi(t)-DikBT∂E∂xiΔt+2DiΔtR+N

For each molecular system, we terminated the simulation when the runs generated at least 600 xk263-HIVp associates and the computed average association time was within the standard deviation of <3%. A successful xk263–HIVp association is defined as xk263 reaches within 11.5 Å ASP25, a residue in the active site of HIVp, and the association time is recorded. Notably, after generating 600 successful bindings, the computed average associating time for each system is within 3%.

In the system, xk263 starts diffusing from the yz-plane at the -x-direction boundary (Figure 1A). A run is ended when xk263 reaches the binding pocket or the +x-direction boundary, and another new run will be started by randomly placing xk263 in a new position on the yz-plane (Figure 3). The clock random number seed is used, so no two runs can be identical. In addition, we can save trajectories of our BD runs to visualize ligand diffusion and binding/unbinding events.

Schematic top view of the simulation box. HIVp is placed on CH3-SAM, and two initial starting positions of xk263, A and B, are on the yz-plane. The simulation is terminated after xk263 binding to HIVp, and another new replica will start. If xk263 reaches the +x-direction wall, the simulation is also terminated, and a new replica will start in a different position on the yz-plane.

 Calculations for System Analysis

For the calculations of diffusion flux, diffusion coefficients, CH3-SAM interaction, and xk263 distribution, we used the trajectories with no termination at the +x-direction boundary. So, the trajectories used for the calculations are continuous and 0.3-μs long for each case.

xk263 Diffusion Flux

The diffusion flux is the rate of molecules transferred across the plane per unit time per unit area, representing the flow of xk263. It was calculated by using the continuous BD trajectories and reported in units of molecules/s.m2. The plane is an imaginary plane at the +x- directional boundary perpendicular to the direction of flow. The calculated flux densities according to the N values (in Equation 2) are reported (Supplementary Table 1).

 Diffusion Coefficient

The diffusion coefficient was calculated by Einstein's relation, < r2> = 2nDt, where r is the displacement in time t, n is the dimensionality, and D is the diffusion coefficient.

 CH<sub>3</sub>-SAM–xk263 Interaction Energy

To calculate the interaction energetics of xk263 with CH3-SAM, we used the selected portion of continuous trajectories with xk263 diffusions within 25 Å above the CH3-SAM. The interaction energy includes 12-6 Lennard-Jones potential and screened Coulombic potential. The cutoffs for vdW and electrostatic potential were 12 Å and 40 Å, respectively.

The calculation of 12-6 Lennard-Jones potential is as follows:

Evdw=4∈(σ12rij12-σ6rij6)

where rij is the distance between atoms i and j, ϵ is the pairwise well-depth parameter and σ is the distance at which the atomic radii meet and where the potential changes sign.

σ=σi+σj2 ϵ=ϵiϵj

The calculation of screened Coulombic potential is as follows:

ESC=keqiqje-krijrijϵ

where rij is the distance between two point-charges i and j from the ligand and receptor, respectively, qi and qj are the partial electron charges of the ligand point charges i and receptor point charges j, k is the screening parameter, ke is the Coulomb constant, kb is the Boltzmann constant and ϵ is the solution dielectric constant.

 xk263 Distribution

The distribution of xk263 in the system was calculated by using continuous BD trajectories. The volume of the system was divided into 11 slices of 20 Å height each. The total number of appearances of xk263 replicas was computed for each slice in 0.3 μs. The calculation of g(z), distribution, is according to Equation 8

g(z)=NM=NnR×vV×f

where N is the total number of appearances of xk263 replicas in a slice within a certain number of frames, M is the expected total number of appearances approximated from a random diffusional motion, nR is the number of xk263 replicas in the system, v is the volume of the slice, V is the total volume of the system, and f is the number of frames to be analyzed.

 Results and Discussion

Understanding ligand–receptor associations in various environments is of vital importance in drug binding in cells or in different experiments. In this study, we analyzed the effects of the SAM surface, the diffusion flux of ligand xk263 and the position of HIVp from the surface and in comparison with results with a static environment (Table 1). HIVp has flexible flaps, which are mostly in closed conformations (Katoh et al., 2003; Kang et al., 2011). To bind with a ligand, the flaps may open spontaneously or are induced to open by the ligand. Because this study focuses on the initial molecular encounter processes, we chose a flap open conformation and used rigid-body BD simulations to model the ligand association. In addition, the orientation of HIVp with respect to the +x direction diffusion flux of xk263 can affect ligand association. Therefore, we examined two HIVp orientations, perpendicular and parallel with the ligand flux (Figure 2). We modeled 60 different systems (Figure 4). We chose the values of xk263 diffusion flux that gradually show the changes of the modeled values with the flux.

Summary of molecular systems used in this study.

Details of model system
HIVp orientations 2
Periodic wall y
Flow direction x
Top Closed in z-direction
Ligand concentration 110 μM
Ligand initial position As a yz-plane
Surface with and without CH3-SAM
Protein height 0 and 20 Å above the surface
Box dimension 400 × 400 × 220 Å3
Termination Reaching +x-direction wall or within 11.5 Å of ASP 25′

Sixty different systems used in this study. The two HIVp orientations are shown in Figure 2.

 Environmental Factors Involved in Ligand Diffusion and Distribution Intermolecular Interactions Between xk263 and CH<sub>3</sub>-SAM

We first examined the interactions between xk263 and CH3-SAM when the ligand diffusion flux varies. Table 2 shows weak van der Waals (vdW) interactions between xk263 and CH3-SAM regardless of the flux, which resulted in no significant adsorption of xk263 on SAM (Figure 6A). Because of the hydrophobic nature of CH3-SAM and xk263, it is not surprising that the electrostatic interaction was close to zero.

Interaction energies between xk263 and CH3-SAM in 7 different xk263 diffusion fluxes.

xk263 diffusion flux (× 1022/s.m2) Eelec (kcal/mol) EvdW (kcal/mol)
0 −0.045 ± 0.14 −1.126 ± 0.89
0.45 −0.056 ± 0.12 −1.098 ± 0.65
1.09 −0.010 ± 0.11 −1.547 ± 1.14
3.43 0.015 ± 0.14 −1.194 ± 0.95
5.62 0.003 ± 0.11 −1.488 ± 0.96
11.29 0.006 ± 0.10 −1.273 ± 1.13
17.20 0.005 ± 0.13 −1.069 ± 0.83

Notably, our BD movement yielded the translational diffusion coefficient of xk263, Dlig = 5.58 ± 0.33 × 10−6 cm2/s when a system had no SAM and no ligand diffusion flux. The modeled ligand diffusion coefficient is in good agreement with analytical values, Danalytical = 4.91 × 10−6 cm2/s, which also validates our simulation setting and BD algorithm (Supplementary Methods). Although the intermolecular vdW was weak, the attractions still affected the diffusion of xk263 molecules, and the molecule diffused a little slower in the presence of SAM with or without xk263 diffusion flux (Figure 5, Supplementary Table 1). Because the flow rate is inversely proportional to the viscosity (Pfitzner, 1976), the diffusion coefficient of xk263 increased with increasing xk263 diffusion flux, as anticipated. However, because the water flow rate was not equal to the xk263 diffusion flux, we did not see a simple linear relationship between the two values.

Plot of translational diffusion coefficient of xk263. Diffusion coefficient increases with xk263 diffusion flux. (see Supplementary Table 1 for numerical data).

 Distribution of xk263 in the Simulation Tube

We further investigated whether the environment may affect the concentration gradient of xk263 along the z-direction when the molecule diffuses in the square tube. The system was partitioned into slices of 20 Å each along z-direction, and Figure 6 shows the distribution of xk263 as a function of their distance above the surface. The distribution of xk263 for a slice, g(z), is the ratio of number of observed xk263 and number of expected xk263 approximated from a random diffusional motion. The systems include HIVp with one orientation shown in Figure 2A. When there was no SAM and the xk263 diffusion flux was small (green line in Figure 6A), the distribution of xk263 was the same throughout the z-direction (distribution ratio = ~1). In contrast, when CH3-SAM was present (Figures 6A–C), xk263 was double that in the first slide (within 20 Å of the SAM) when diffusion flux of xk263 was <1.09 × 1022 molecules/s m2. The results suggest that even if the intermolecular attraction between xk263 and CH3-SAM is weak, the concentration can be increased by 2-fold, which helps to increase the probability of a molecular encounter. Studies for various systems also showed that local molecular concentration can be effected by a surface; for example, surface catalyzed biomolecular reaction using nanoreactors, nucleotidase co-localization, and absorption of biomacromolecules on hydrogels (Roa et al., 2017; Pérez-Mas et al., 2018; Rahmaninejad et al., 2020).

Plot of concentration distribution of xk263, g(z), as a function of the distance above the CH3-SAM, z, in seven different diffusion fluxes (A–G). The space was partitioned into equally sized slices with a height of 20 Å. g(z) is the ratio of number of observed xk263 and the number of expected xk263 approximated from a random diffusional motion.

The SAM retains xk263 near the surface only when the diffusion flux of the ligand was small. When the flux increases, the local concentration of xk263 near the surface decreased, and xk263 had the highest concentration near the middle of the z-direction (see ~100 Å in Figures 6D–G). The distribution of concentration gradient along the z-direction was the same with or without SAM, which suggests that the weak intermolecular interactions can be altered easily when xk263 has diffusion flux in the x-direction. Our results are consistent with experiments showing the quadratic curve of concentration gradients with flow (Chen et al., 2017). With larger diffusion flux, a significantly smaller amount of xk263 diffused near the surface where HIVp locates, and the reduced distribution also contributed toward probability of small xk263–HIVp association, discussed later.

 HIVp–xk263 Association in Different Environments HIVp Immobilized on the Surface

We modeled the average xk263–HIVp association time and the binding percentage when HIVp was 0 Å above the surface with SAM (Tables 3, 4 with SAM) or without SAM. Theoretically, if xk263 employs only 3-dimensional diffusion within a sphere, the average association time to bind to HIVp located in the center of the sphere is 2648.4 ns when the system has no external flux (Supplementary Methods) (Adam and Delbrück, 1968). The association time modeled by our simulation is 1491.55 ns is close to the theoretical value. The faster binding time than theory is due to use of a box, and xk263 did not need to search the whole sphere to associate with HIVp. Notably, existing association rate theories describe system without diffusion flux (Berg and Purcell, 1977; Szabo et al., 1980; Shoup and Szabo, 1982) and further theoretical work in molecular association with buffer flowing is needed. Although the interactions between xk263 and CH3-SAM were small, to examine the environment effects, we turned off the intermolecular vdW and electrostatic interactions, so the system was in a confined space without intermolecular attractions or repulsion (without SAM). The x-direction diffusion flux of xk263 gradually increased from zero to 1.7 × 1023 molecules/s m2, and the average xk263–HIVp association time decreased non-linearly. The correlation was observed in a previous study, showing that the association time for human immunoglobulin G (IgG) and Staphylococcus aureus protein A (SpA) decreased with increasing flow rate (Ogi et al., 2008). The average association time did not differ without or with SAM in the presence of diffusion flux (Table 3). This is not surprising because the CH3-SAM provides only weak attractions to xk263, and the surface neither increased the local concentration of xk263 to assist ligand association nor retained the ligand from binding to HIVp (Figure 6) (Cholko et al., 2020). However, if the attraction between a ligand and surface is large, it may affect experimental results. As a result, existing experimental studies considered the choice of the surface for the sensor or SPR to optimize experimental sensitivity and accuracy. For example, phospholipid/alkane bilayers might be a better option for molecular binding in biomolecular systems using SPR (Plant et al., 1995). These phospholipid derivatized surfaces may bring non-specific interactions for xk263 which is absent in our current model, which can result in slightly longer association time.

Average association time (tavg) with different xk263 diffusion fluxes.

HIVp perpendicular to the yz–plane, tavg (ns) HIVp parallel to the yz-plane, tavg (ns)
xk263 diffusion flux (× 1022 /s.m2) 0Å above the surface, With SAM 0Å above the surface, No SAM 20Å above the surface, With SAM 20Å above the surface, No SAM 0Å above the surface, With SAM 0Å above the surface, No SAM 20Å above the surface, With SAM 20Å above the surface, No SAM
0 1491.55 1064.41 1213.45 901.99 1405.77 1084.49 1124.5 905.21
0.45 126.50 122.00 108.70 110.17 145.54 124.90 128.30 111.68
1.09 55.56 51.93 43.61 45.28 53.71 50.05 47.98 44.65
3.43 18.73 18.83 15.042 15.57 17.77 17.83 16.23 15.53
5.62 11.08 10.48 10.31 9.47 10.78 11.01 9.75 9.55
11.29 5.62 5.52 5.52 4.80 5.53 5.38 4.89 4.89
17.20 3.69 3.59 3.33 3.25 3.58 3.69 3.27 3.26

Inverse of average association time per molarity and percentage of successful binding with different xk263 diffusion fluxes.

xk263 diffusion flux (×1022 /s.m2) 0Å above the surface, With SAM 0Å above the surface, No SAM 20Å above the surface, With SAM 20Å above the surface, No SAM
1/tavg perM (/nsM) Association percentage 1/tavg per M (/nsM) Association percentage 1/tavg per M (/nsM) Association percentage 1/tavg per M (/nsM) Association percentage
HIVp perpendicular to the yz-plane
0 6.09 18.59 8.54 18.14 7.49 18.67 10.08 19.18
0.45 71.86 3.02 74.52 3.82 83.63 3.53 82.52 4.24
1.09 163.63 1.22 175.07 1.73 208.45 1.67 200.75 2.16
3.43 485.44 0.52 482.66 0.64 604.35 1.00 583.76 0.89
5.62 820.41 0.35 867.45 0.33 881.41 0.90 959.31 0.65
11.29 1618.17 0.18 1646.61 0.14 1646.43 0.70 1892.99 0.43
17.20 2461.66 0.11 2532.29 0.07 2732.55 0.46 2793.33 0.34
HIVp parallel to the yz-plane
0 6.47 17.15 8.38 17.85 8.08 18.99 10.04 20.61
0.45 62.46 3.25 72.78 3.88 70.86 3.43 81.40 4.22
1.09 169.25 1.21 181.64 1.71 189.48 1.86 203.61 1.92
3.43 511.39 0.38 509.78 0.48 559.99 0.71 585.19 0.75
5.62 842.87 0.22 825.62 0.21 931.88 0.40 952.29 0.42
11.29 1645.34 0.05 1690.36 0.05 1858.74 0.14 1858.13 0.15
17.20 2538.08 0.02 2461.19 0.02 2781.88 0.06 2788.36 0.07

We report the inverse of association time per molarity and the association percentage (Table 4). When xk263 was freely diffused in the tube without diffusion flux, the association rate approximated using the average association time, 1/1491.55 ns per M, yielded kon = 7.8 × 109 1/Ms. Using the diffusion coefficient from the BD run (Supplementary Table 2), our simulated kon was approximately a half of the analytical value, kon = 4πRD = 1.2 × 1010 1/Ms because the analytical formula does not consider molecular geometry in association. Our modeled kon was also slightly slower than that measured using SPR for another highly similar cyclic urea compound, 2.5 × 1010 1/Ms (Markgren et al., 2002). The first phase in SPR utilizes buffer flowing and detects binding signals, which require successful collision to have the correct ligand and HIVp orientation and enough energy. The flow may influence successful collision, especially on biomolecular systems which usually have complex molecular geometry. Of note, the inverse association time increased linearly with xk263 diffusion flux (Supplementary Figure 1), which suggests that the association rate may also increase linearly with the flow. However, our modeled diffusion coefficients (Supplementary Table 2) show an exponential increase with the xk263 diffusion flux (Figure 5). When particles have flux moving along the +x direction, the ligand did not freely diffuse in all directions, and the equation for the diffusion-controlled association system, kon = 4πRD, was no longer suitable to approximate the association rate constants. Therefore, we cannot directly estimate kon from the increased diffusion coefficient.

The additional +x direction drift velocity also reduced the initial encounter probability. When the flux was >3.43 × 1022 molecules/s.m2, xk263 passed the tube without sufficient time to freely diffuse in the y and z directions, which resulted in reduced search space and tremendously decreased association percentage. Large diffusion flux also yielded fewer xk263 molecules staying near both the SAM and HIVp surface. Notably, xk263 may diffuse on the surface of HIVp before binding to the pocket. The xk263 diffusion flux in the +x direction perturbed the interactions between xk263 and HIVp, and the weakened intermolecular attractions also led to decreased association percentage.

In typical experimental settings, an immobilized protein can have multiple orientations. Therefore, we chose two representative orientations—the binding pocket of HIVp perpendicular or parallel to the flux of xk263—to examine whether the orientation against the flow affects ligand binding. Tables 3, 4, Figures 7, 8 show no difference in average ligand association time and binding percentage with the two orientations. Because the binding pocket of HIVp is widely accessible for the ligand, the impact of the orientation was insignificant. However, some proteins exhibit a certain direction for ligands to enter the binding pocket, and their orientation against the flux of the ligands can affect the ligand–protein association.

Plot of average association time vs. xk263 diffusion flux. HIVp is placed (A) perpendicular to the yz-plane and (B) parallel to the yz-plane; Similar data imply no effect of HIVp orientation on the association time. In each plot, coinciding lines (red, green, blue, and magenta) show no effect of hydrophobic CH3-SAM or height of the receptor from the surface.

Plot of association percentage vs. xk263 diffusion flux under two different HIVP's orientation (A,B). Association percentage is defined as the fraction of successful binding from the total number of trajectories in the system.

 Artificially Localized HIVp 20 Å Above the Surface

Although different HIVp orientations did not affect xk263 association time, existing studies showed that the dimer interface (the bottom part of HIVp) is the most popular region for ligands' initial encounter of HIVp (Roberts and Chang, 2016). After reaching the dimer interface, the ligand can undergo surface diffusion to reach the binding pocket. When we placed HIVp 0 Å above the surface, this bottom region became partially inaccessible to xk263. Therefore, immobilized HIVp was artificially placed 20 Å above the surface so that xk263 can easily reach the bottom region. When the whole HIVp surface was accessible to xk263, the association time was reduced by 10–20% when the xk263 diffusion flux was <1.09 × 1022 molecules/s.m2, which suggests that the additionally available bottom region accelerates xk263 binding (Table 3) and slightly increases the association percentage (Table 4). Even when HIVp was placed 20 Å above the surface, the location was not far enough to eliminate the intermolecular xk263–SAM interactions. As a result, xk263 spent longer time near SAM when xk263 diffusion flux was small, thus resulting in longer time to associate with HIVp when SAM was present rather than absent. When the flux exists, the association time was the same when the SAM was present or absent, regardless the position of HIVp (Table 3).

 Conclusion

Ligand–protein encounters can occur in any environment that may provide additional intermolecular interactions or the transporting forces to the ligand. For example, in experimental settings such as using SPR to study ligand–protein binding, the choice of the surface and the flow rate are all optimized for measurements. In this study, we used BD simulations to investigate the effect of the CH3-SAM surface on xk263–HIVp association. The non-polar surface provided weak intermolecular vdW attractions with xk263 to slightly increase the ligand binding time. The effects quickly vanished when xk263 had an x-direction diffusion flux, and the association time decreased when the diffusion flux increased regardless of the presence of SAM or only a special plane. With no diffusion flux, xk263 was twice more concentrated within 20 Å of the SAM surface because of the xk263–SAM interactions. The concentration gradient did not increase binding time but increased xk263–HIVp binding probability. When the xk263 diffusion flux increased, the middle region of the square tube had the highest xk263 concentration, which is the same as existing experiments showing a quadratic curve of concentration gradients with flow. The results also show that when HIVp was placed on the surface (0 Å above the surface), the bottom part of the protein, a known high-probability site of the xk263–HIVp first encounter, was not accessible to xk263. Because xk263 could bind non-specifically on the HIVp surface and then utilize surface diffusion to reach the binding site, occluding this bottom region increased ligand binding time by 10–20% in the static environment. However, with large xk263 diffusion flux, all the weak intermolecular interactions and searching along the HIVp surface were eliminated, which resulted in a fast association time and small binding probability. This work brings insights into how ligand diffusion flux and the environment may affect ligand–protein association.

 Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

 Author Contributions

SK ran and analyzed the simulations, produced figures, and wrote the manuscript. C-eC designed experiments, analyzed the simulation, and wrote the manuscript. Both authors contributed to the article and approved the submission version.

 Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

We thank Dr. Jason Quan Cheng for helpful discussions regarding SPR experiments. We thank Tim Cholko for help in modifying the GeomBD program and Jianan Sun for providing HIVp open structures.

 Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2021.659687/full#supplementary-material

 References Adam G. Delbrück M. (1968). Reduction of dimensionality in biological diffusion processes. Struct. Chem. Mol. Biol. 198, 198–215. Andrews S. S. (2017). Smoldyn: particle-based simulation with rule-based modeling, improved molecular interaction and a library interface. Bioinformatics 33, 710–717. 10.1093/bioinformatics/btw70028365760 Baron R. McCammon J. A. (2013). Molecular recognition and ligand association. Annu. Rev. Phys. Chem. 64, 151–175 10.1146/annurev-physchem-040412-110047 Berg H. C. Purcell E. M. (1977). Physics of chemoreception. Biophys. J. 20, 193–219. 10.1016/S0006-3495(77)85544-6 Bernetti M. Masetti M. Rocchia W. Cavalli A. (2019). Kinetics of drug binding and residence time. Annu. Rev. Phys. Chem. 70, 143–171. 10.1146/annurev-physchem-042018-05234030786217 Britton J. Majumdar S. Weiss G. A. (2018). Continuous flow biocatalysis. Chem. Soc. Rev. 47, 5891–5918. 10.1039/C7CS00906B Chang C.-E. Shen T. Trylska J. Tozzini V. McCammon J. A. (2006). Gated binding of ligands to HIV-1 protease: brownian dynamics simulations in a coarse-grained model. Biophys. J. 90, 3880–3885. 10.1529/biophysj.105.07457516533835 Chen X. Hu Z. Zhang L. Yao Z. Chen X. Zheng Y. . (2017). Numerical and experimental study on a microfluidic concentration gradient generator for arbitrary approximate linear and quadratic concentration curve output. Int. J. Chem. Reactor Eng. 16. 10.1515/ijcre-2016-0204 Cholko T. Barnum J. Chang C. E. (2020). Amyloid-β (Aβ42) peptide aggregation rate and mechanism on surfaces with widely varied properties: insights from brownian dynamics simulations. J. Phys. Chem. B 124, 5549–5558. 10.1021/acs.jpcb.0c0292632525673 Cholko T. Kaushik S. Chia-en A. C. (2019). Dynamics and molecular interactions of single-stranded DNA in nucleic acid biosensors with varied surface properties. Phys. Chem. Chem. Phys. 21, 16367–16380. 10.1039/C9CP02441G31309941 De Biase P. M. Markosyan S. Noskov S. (2015). BROMOC suite: Monte Carlo/Brownian dynamics suite for studies of ion permeation and DNA transport in biological and artificial pores with effective potentials. J. Comput. Chem. 36, 264–271. 10.1002/jcc.2379925503688 Di Cera E. (2017). Mechanisms of ligand binding. Biophys. Rev. 1:011303. 10.1063/5.0020997 Dill K. Bromberg S. (2012). Molecular driving Forces: Statistical Thermodynamics in Biology, Chemistry, Physics, and Nanoscience. New York, NY: Garland Science. 10.4324/9780203809075 Długosz M. Zieliński P. Trylska J. (2011). Brownian dynamics simulations on CPU and GPU with BD_BOX. J. Comput. Chem. 32, 2734–2744. 10.1002/jcc.2184721638295 Ershov P. V. Mezentsev Y. V. Kaluzhskiy L. A. Ivanov A. S. (2020). Phenanthridine derivatives as potential HIV-1 protease inhibitors. Biomed. Rep. 13:66. 10.3892/br.2020.137333149910 Ghosh A. K. Osswald H. L. Prato G. (2016). Recent progress in the development of HIV-1 protease inhibitors for the treatment of HIV/AIDS. J. Med. Chem. 59, 5172–5208. 10.1021/acs.jmedchem.5b0169726799988 Hinman S. S. McKeating K. S. Cheng Q. (2018). Surface plasmon resonance: material and interface design for universal accessibility. Anal. Chem. 90:19. 10.1021/acs.analchem.7b0425129053253 Huang X. Britto M. D. Kear-Scott J. L. Boone C. D. Rocca J. R. Simmerling C. . (2014). The role of select subtype polymorphisms on HIV-1 protease conformational sampling and dynamics. J. Biol. Chem. 289, 17203–17214. 10.1074/jbc.M114.57183624742668 Huang Y. M. Raymundo M. A. Chen W. Chang C. A. (2017). Mechanism of the association pathways for a pair of fast and slow binding ligands of HIV-1 protease. Biochemistry 56, 1311–1323. 10.1021/acs.biochem.6b0111228060481 Huber G. A. McCammon J. A. (2010). Browndye: a software package for Brownian dynamics. Comput. Phys. Commun. 181, 1896–1905. 10.1016/j.cpc.2010.07.02221132109 Huber G. A. McCammon J. A. (2019). Brownian dynamics simulations of biological molecules. Trends Chem. 1, 727–738. 10.1016/j.trechm.2019.07.00832309795 Kang M. Roberts C. Cheng Y. Chang C. E. (2011). Gating and intermolecular interactions in ligand-protein association: coarse-grained modeling of HIV-1 protease. J. Chem. Theory Comput. 7, 3438–3446. 10.1021/ct200488526598172 Katoh E. Louis J. M. Yamazaki T. Gronenborn A. M. Torchia D. A. Ishima R. (2003). A solution NMR study of the binding kinetics and the internal dynamics of an HIV-1 protease-substrate complex. Protein Sci. 12, 1376–1385. 10.1110/ps.030070312824484 Kohl N. E. Emini E. A. Schleif W. A. Davis L. J. Heimbach J. C. Dixon R. . (1988). Active human immunodeficiency virus protease is required for viral infectivity. Proce. Natl Acad. Sci. U.S.A. 85, 4686–4690. 10.1073/pnas.85.13.46863290901 Lam P. Jadhav P. K. Eyermann C. J. Hodge C. N. Ru Y. Bacheler L. T. . (1994). Rational design of potent, bioavailable, nonpeptide cyclic ureas as HIV protease inhibitors. Science 263, 380–384. 10.1126/science.82788128278812 Li D. Liu M. S. Ji B. Hwang K. C. Huang Y. (2012). Identifying the molecular mechanics and binding dynamics characteristics of potent inhibitors to HIV-1 protease. Chem. Biol. Drug Design 80, 440–454. 10.1111/j.1747-0285.2012.01417.x22621379 Lin Y.-C. Kim W. K. Dzubiella J. (2020). Coverage fluctuations and correlations in nanoparticle-catalyzed diffusion-influenced bimolecular reactions. J. Phys. Chem. C 124, 24204–24214. 10.1021/acs.jpcc.0c06898 Madura J. D. Briggs J. M. Wade R. C. Davis M. E. Luty B. A. Ilin A. . (1995). Electrostatics and diffusion of molecules in solution: simulations with the University of Houston Brownian Dynamics program. Comput. Phys. Commun. 91, 57–95. 10.1016/0010-4655(95)00043-F Markgren P.-O. Schaal W. Hämäläinen M. Karlén A. Hallberg A. Samuelsson B. . (2002). Relationships between structure and interaction kinetics for HIV-1 protease inhibitors. J. Med. Chem. 45, 5430–5439. 10.1021/jm020837012459011 Martinez M. Bruce N. J. Romanowska J. Kokh D. B. Ozboyaci M. Yu X. . (2015). SDA 7: A modular and parallel implementation of the simulation of diffusional association software. J. Comput. Chem. 36, 1631–1645. 10.1002/jcc.2397126123630 McCammon J. A. (2011). Gated diffusion-controlled reactions. BMC Biophys. 4, 1–5. 10.1186/2046-1682-4-4 Northrup S. Laughner T. Stevenson G. (1999). MacroDox Macromolecular Simulation Program. Cookeville, TN: Tennessee Technological University, Department of Chemistry. Northrup S. H. Allison S. A. McCammon J. A. (1984). Brownian dynamics simulation of diffusion-influenced bimolecular reactions. J. Chem. Phys. 80, 1517–1524. 10.1063/1.446900 Ogi H. Fukunishi Y. Omori T. Hatanaka K. Hirao M. Nishiyama M. (2008). Effects of flow rate on sensitivity and affinity in flow injection biosensor systems studied by 55-MHz wireless quartz crystal microbalance. Anal. Chem. 80, 5494–5500. 10.1021/ac800459g18553940 Ozbabacan S. A. Gursoy A. Keskin O. Nussinov R. (2010). Conformational ensembles, signal transduction and residue hot spots: application to drug discovery. Curr. Opin. Drug Discov. Devel. 13, 527–537.20812144 Pang X. Zhou H.-X. (2017). Rate constants and mechanisms of protein–ligand binding. Ann. Rev. Biophys. 46, 105–130. 10.1146/annurev-biophys-070816-03363928375732 Pérez-Mas L. Martín-Molina A. Quesada-Pérez M. Moncho-Jordá A. (2018). Maximizing the absorption of small cosolutes inside neutral hydrogels: steric exclusion versus hydrophobic adhesion. Phys. Chem. Chem. Phys. 20, 2814–2825. 10.1039/C7CP07679G29323684 Pfitzner J. (1976). Poiseuille and his law. Anaesthesia 31, 273–275. 10.1111/j.1365-2044.1976.tb11804.x779509 Planchestainer M. Contente M. L. Cassidy J. Molinari F. Tamborini L. Paradisi F. (2017). Continuous flow biocatalysis: production and in-line purification of amines by immobilised transaminase from Halomonas elongata. Green Chem. 19, 372–375. 10.1039/C6GC01780K Plant A. L. Brighamburke M. Petrella E. C. Oshannessy D. J. (1995). Phospholipid/alkanethiol bilayers for cell-surface receptor studies by surface plasmon resonance. Anal. Biochem. 226, 342–348. 10.1006/abio.1995.12347793636 Prabowo B. A. Purwidyantri A. Liu K.-C. (2018). Surface plasmon resonance optical sensor: a review on light source technology. Biosensors 8:80. 10.3390/bios803008030149679 Rahmaninejad H. Pace T. Bhatt S. Sun B. Kekenes-Huskey P. (2020). Co-localization and confinement of ecto-nucleotidases modulate extracellular adenosine nucleotide distributions. PLoS Comput. Biol. 16:e1007903. 10.1371/journal.pcbi.100790332584811 Roa R. Kim W. K. Kanduc M Dzubiella J. Angioletti-Uberti S. (2017). Catalyzed bimolecular reactions in responsive nanoreactors. ACS Catal. 7, 5604–5611. 10.1021/acscatal.7b0170128966839 Roberts C. C. Chang C.-E. (2016). Analysis of ligand–receptor association and intermediate transfer rates in multienzyme nanostructures with all-atom brownian dynamics simulations. J. Phys. Chem. B 120, 8518–8531 10.1021/acs.jpcb.6b0223627248669 Roberts C. C. Chang C. E. (2015). Modeling of enhanced catalysis in multienzyme nanostructures: effect of molecular scaffolds, spatial organization, and concentration. J. Chem. Theory Comput. 11, 286–292 10.1021/ct500748226574226 Saadat A. Khomami B. (2015). Matrix-free Brownian dynamics simulation technique for semidilute polymeric solutions. Phys. Rev. E 92:033307 10.1103/PhysRevE.92.03330726465586 Schöneberg J. Noé F. (2013). ReaDDy-a software for particle-based reaction-diffusion dynamics in crowded cellular environments. PLoS ONE 8:e74261(e2013). 10.1371/journal.pone.007426124040218 Shoup D. Szabo A. (1982). Role of diffusion in ligand binding to macromolecules and cell-bound receptors. Biophys. J. 40, 33–39. 10.1016/S0006-3495(82)84455-X7139033 Spinelli S. Liu Q. Alzari P. Hirel P. Poljak R. (1991). The three-dimensional structure of the aspartyl protease from the HIV-1 isolate BRU. Biochimie 73, 1391–1396. 10.1016/0300-9084(91)90169-21799632 Szabo A. Schulten K. Schulten Z. (1980). First passage time approach to diffusion controlled reactions. J. Chem. Phys. 72, 4350–4357. 10.1063/1.439715 Szabo A. Shoup D. Northrup S. H. McCammon J. A. (1982). Stochastically gated diffusion-influenced reactions. J. Chem. Phys. 77, 4484–4493. 10.1063/1.444397 Tozzini V. McCammon J. A. (2005). A coarse grained model for the dynamics of flap opening in HIV-1 protease. Chem. Phys. Lett. 413, 123–128. 10.1016/j.cplett.2005.07.07517029846 Votapka L. W. Jagger B. R. Heyneman A. L. Amaro E. R. (2017). SEEKR: simulation enabled estimation of kinetic rates, a computational tool to estimate molecular kinetics and its application to trypsin–benzamidine binding. J. Phys. Chem. B 121, 3597–3606. 10.1021/acs.jpcb.6b0938828191969 Zotter A. Bäuerle F. Dey D. Kiss V. Schreiber G. (2017). Quantifying enzyme activity in living cells. J. Biol. Chem. 292, 15838–15848. 10.1074/jbc.M117.79211928784664

Funding. This study was supported by the US National Institutes of Health (GM-109045), US National Science Foundation (MCB-1932984), and NSF national supercomputer centers (TG-CHE130009).