Dynamic model of estrogen docking using multiscale analysis

Abstract

This work investigates the application of a new multiscale analysis to dynamic modeling of the interaction between the estrogen molecule and its receptor. The key notion here is to explore the applications of this approach to objects as small as the estrogen molecule, a study that has not yet been done. This work lays the foundation for the development of a new theoretical screening technique to identify carcinogens. In this work, a three-dimensional, coarse grained approximation of estrogen is modeled. In order to facilitate the application of the techniques used herein, several reasonable simplifications of the estrogen model have been made. It has been observed that the time required to numerically integrate the classical Newton–Euler model of estrogen is quite long, because a small time step size, on the order of 0.1 ps ($10^{-13}$ s) must be used to capture the dynamics. Here, a new multiscale analysis is used to develop a scaled model that can be numerically integrated in less time, yet accurately predict the system’s behavior .

Appendices

Appendix 1: Viscous forces

In order to determine the viscous friction properties, it is necessary to consider the characteristics of the fluid surrounding the small particle. The key issue here is determining whether the fluid should be modeled as a continuum, or as discrete, individual molecules; the goal is not to model individual fluid molecules. However, there is difficulty in determining a coefficient of viscous friction for the complex shape of estrogen. Thus, the effort is to apply viscous friction forces to each atom and sum them to find the resultant force on the coarse grained estrogen. This allows determination of a overall coefficient of viscous friction for the system, which enables the multiscale analysis in Sect.4.

The first step in the process is to determine whether the continuum modeling approach can be used with the coarse grained estrogen model. This can be investigated using the Knudsen number, $Kn$. It is the ratio of the fluid mean free path and the characteristic length of system:

$$\begin{aligned} \displaystyle Kn \ = \ \frac{\uplambda _\mathrm{water}}{L_\mathrm{estrogen}} \ = \ \frac{ 0.25\,\hbox {nm}}{1.18\,\hbox {nm}} \ = \ 0.21 \end{aligned}$$

(17)

where $\uplambda _\mathrm{water}$ is the mean free path of water and $L_\mathrm{estrogen}$ is the characteristic length of estrogen. Since $Kn<1$, the viscous friction force can be modeled using the continuum formulation. Thus, it should be reasonable to use Stoke’s Law to calculate drag coefficients.

However, because the estrogen is small, it is unclear how its surface interacts with the surrounding fluid to create friction or drag. This is referred to as the no-slip boundary condition, which indicates whether the fluid sticks to the estrogen’s surface creating larger drag forces, or slipping occurs between the fluid and the estrogen, creating less drag. This has been verified experimentally and theoretically [24]. This condition can be checked by calculating a slip correction factor, $C_c$, [3, 13],

$$\begin{aligned} \displaystyle C_c = 1 + 2 Kn \left( 1.257 + 0.4e^{\frac{-1.1}{2Kn}} \right) = 1.6 \end{aligned}$$

(18)

which divides the drag coefficients.

Since it is difficult to determine the drag on the complex estrogen shape, here an approximation is used by portioning the resultant viscous drag force to each atom comprising estrogen. Assuming each atom is spherical makes it simple to calculate the viscous drag coefficient. This is then multiplied by the correction factor $C_c$ to obtain the coefficient of viscous friction for an atom. The drag force on each atom is included in the equations of motion where they sum to create the generalized active force described by $\varvec{\varGamma }_F = \beta D \dot{\mathbf{q}}$. Application of the drag forces to each atom is illustrated in Fig. 7a, which can be replaced by an equivalent force system consisting of the resultant drag force and moment at the mass center, as shown in Fig. 7b.

Fig. 7

Drag forces and moments: a on individual atoms, b resultant drag forces and moments

The linear, $\beta _v$, and rotational, $\beta _w$, drag coefficients for a sphere are given by

$$\begin{aligned} \bar{\beta }_v = 6 \pi \mu r,\quad \bar{\beta }_{w} = 8 \pi \mu r^3 \end{aligned}$$

(19)

where $\mu$ is the Viscosity of water and $r$ is the radius of the sphere. Here, the fluid is assumed to be a Newtonian fluid with uniform viscosity at $20\,^{\circ }$C. The modified drag coefficients are found as:

$$\begin{aligned} \beta _v = \frac{\bar{\beta }_v}{C_c},\quad \beta _w = \frac{\bar{\beta }_w}{C_c} \end{aligned}$$

(20)

The viscous friction force and moment acting on each atom can be calculated from modified Stroke’s law:

$$\begin{aligned} \mathbf{F}_{f_i} = -\beta _v ^{N}\mathbf{V}_i,\quad \mathbf{M}_{f_j} = -\beta _{w} ^{N}{\varvec{\omega }}^j \end{aligned}$$

(21)

where $^{N}\mathbf{V}_i$ and $^{N}{\varvec{\omega }}^j$, $i$ is a point and $j$ is a body, are the translational and rotational velocities of the sphere observed from the inertial reference frame, $N$, in Fig. 2.

In order to find the viscous friction coefficients for the base structure for the coarse grained estrogen model, A1–A19, the resultant friction force and moment can be examined. This is accomplished as follows:

$$\begin{aligned} \displaystyle \mathbf{F}_{f_{A_o}}&= \sum ^{19}_{i=1}{\mathbf{F}_{f_{Ai}}}\end{aligned}$$

(22)

$$\begin{aligned} \mathbf{M}_{f_A}&= \displaystyle \sum ^{19}_{j=1}{\mathbf{M}_{f_{Aj}} \ + \ \mathbf{P}_{A_oAj} \times \mathbf{F}_{f_{Ai}}}\nonumber \\&= \displaystyle \sum ^{19}_{j=1}{-\beta _w ^N{\varvec{\omega }}^A \ + \ \mathbf{P}_{A_oAj} \times \mathbf{F}_{f_{Ai}}} \end{aligned}$$

(23)

where $\mathbf{P}_{ij}$ is the position vector from point $i$ to point $j$. These resultant forces and moments are applied to the coarse grained estrogen model. The unscaled Newton–Euler model is then simulated, and the resulting motion can be used to determine the overall coefficient viscous friction coefficient $\beta$ in (3). Here, the translational drag coefficient was used to estimate the overall $\beta$

$$\begin{aligned} \beta = \frac{ \Vert \mathbf{F}_{f_{A_o}} \Vert }{ 0.001 + \Vert ^{N}\mathbf{V}_{A_o} \Vert } \end{aligned}$$

(24)

The $0.001$ is added to the denominator to prevent an infinity value when the velocity is equal to zero. The mean value of $\beta$ calculated from (24) is 98.52 zg/hs as shown in Fig. 8. The assemblage of these viscous friction forces comprises $\varvec{\varGamma }_{F}$ in the equations of motions.

Fig. 8

$\beta$ Estimate from simulation without Brownian motion

Appendix 2: Charge forces

The docked estrogen forms four Hydrogen bonds with the ER. The following assumptions are made in modeling the charge forces in this model:

• 1. Since the potential of the Hydrogen bonding is complex [21], for simplicity, the Coulomb potential is used.

• 2. Modified hard sphere repulsion is used to model the repulsive forces that exist when the atoms get closer than their Van der Waals radii.

The interaction potential of the Coulomb charge or Coulomb potential, $V_a$, is given by the relation, (25) [21].

$$\begin{aligned} V_a(r) = \frac{Q_1 \ Q_2}{4 \ \pi \ \varepsilon \ r} \end{aligned}$$

(25)

where $V_a(r)$ is the Coulomb Potential, $Q_1$, $Q_2$ are Coulomb Charges, $\varepsilon$ is the permittivity of the medium; $\varepsilon = \varepsilon _o \varepsilon _r$, where $\varepsilon _o$ is permittivity of vacuum and $\varepsilon _r$ is the relative permittivity of medium and $r$ is the distance between charges.

These charges attract when their separation, $r$, is greater than the Van der Waalls radius, and repulsive when less. The repulsive force can be modeled as [21],

$$\begin{aligned} V_r(r) = \frac{B_{i}}{r^6} \end{aligned}$$

(26)

where $V_r$ is the repulsive potential, $B_{i}$ is a constant dependent on Van der Waals radii of the atoms and bond length, and $r$ is the distance between the charges.

The interaction potential, $V$, of a bond is the sum of the Coulomb (25) and repulsive (26) potentials:

$$\begin{aligned} V(r) = V_a(r) + V_r(r) = \frac{e^2}{4 \pi \varepsilon } \left( {\frac{Q_1 Q_2}{r}+\frac{B}{r^6}} \right) \end{aligned}$$

(27)

where $B = B_i \left( 4 \pi \varepsilon /e^{2} \right)$, $e$ is value of unit charge. The force associated with this potential, $F(r)$, can be calculated by taking the derivative of interaction energy with respect to the distance as shown in (28)

$$\begin{aligned} F(r) = \frac{d V(r)}{dr} \end{aligned}$$

(28)

The charge values $Q_1, \ldots , Q_7$ in Fig. 9 are obtained from charge density studies [28] and from the Ab initio quantum mechanical studies [16]. The values used are given in Table 3.

Fig. 9

$\mathbf{F}_1$, $\mathbf{F}_2$, and $\mathbf{F}_3$ represent the force due to Hydrogen bonds and are modeled using (28). The values of charges $Q_1, \ldots , Q_6$ are given in Table 3

The bond distances used in this model are given in Table 3 [28]. The constant $B$ in the repulsive potential in (26) defines the point at which the interaction potential for a bond is minimum. In other words, the force between the atoms is zero at this point, and hence the system is at equilibrium. The value of the constant $B$ is calculated by finding the local minimum of the total interaction potential in (27), when $r$ is equal to the bond distance of the Hydrogen bond. The values used in this study are listed in Table 3.

The force magnitudes for each of the modeled charges are shown in Fig. 10. Examples of the vector forces modeled for points D3 and $B_o$ in Fig. 9 are

$$\begin{aligned} \mathbf{F}_{1_{D3}} \ = \ F_1 \ \frac{ \mathbf{P}_{D3B_o} }{ \Vert \mathbf{P}_{D3B_o} \Vert } \ = \ F( \Vert \mathbf{P}_{D3B_o} \Vert ) \ \frac{ \mathbf{P}_{D3B_o} }{ \Vert \mathbf{P}_{D3B_o} \Vert } \end{aligned}$$

(29)

and the charge force acting on $B_o$ is

$$\begin{aligned} \mathbf{F}_{1_{B_o}} \ = \ -\mathbf{F}_{1_{D3}} \end{aligned}$$

(30)

The collection of these charge forces defines $\varvec{\varGamma }_C$ in the equations of motion.

Fig. 10

Magnitude of charge forces versus the distance between the charges

Table 3 Definition of quantities used

Appendix 3: Brownian motion

Here, the idea is that the Brownian force is expected to have a resultant force and moment as shown in Fig. 7b, and this resultant force can be modeled as a series of forces and moments on individual atoms as shown in Fig. 7a.

The random forces and moments in the model that produce Brownian motion are implemented as Gaussian white noise. They are applied to each atom, similar to Fig. 7a. For example, the random force and moment on a Hydrogen atom is defined as

$$\begin{aligned}&\mathbf{F}_{A1} = B_{A1_{1}}(t) \ \widehat{\mathbf{N}}_1 + B_{A1_{2}}(t) \ \widehat{\mathbf{N}}_2 \nonumber \\&\mathbf{M}_{A1} = \bar{r}_O \ B_{A1_{3}}(t) \end{aligned}$$

(31)

where the $\bar{r}_O$ is a characteristic radius of Oxygen atom. The $B_{oi}(t)$ represents forces produced by randomly fluctuating thermal noise on the atom A1, an Oxygen atom. Each component of the random force and moment is treated independently as a normally distributed random variable. They have the following expectations, E[$\cdots$], or weighted average values,

$$\begin{aligned} E[B_{A1_{i}}(t)] = {\left\langle \ B_{A1_{i}}(t) \ \right\rangle } = 0 = \mu \end{aligned}$$

(32)

and are governed by a fluctuation–dissipation relaxation expressed as

$$\begin{aligned} E \left[ \ B_{A1_{i}}(t_1) \ B_{A1_{j}}(t_2)\ \right] = 2 \ (\beta _v)_O \ k_B \ \mathrm{T} \ \delta ({t}_{i}-{t}_{j}) \delta _{i,j} \end{aligned}$$

(33)

where $k_B$ and T are the Boltzmann constant and absolute temperature [20]. The relations in (33) imply that there is no time dependency between the random process over time; the random sequence of forces does not repeat regularly. In addition, (32) and (33) imply

$$\begin{aligned} E[ {B}_{A1_{i}}^{2}(t) ] = 2 \ \beta \ k_B \ \mathrm{T} = \mathrm{Var}(B_{A1_{i}}(t)) = {\sigma }^{2} \end{aligned}$$

(34)

which are the variances of the $B_{A1_{i}}$. For Carbon $\sigma = 41.81 \times 10^{-2}$ and for Oxygen $\sigma = 36.05 \times 10^{-2}$ while the mean for both is $\mu =0$. Thus, the $B_{A1_{i}}$ can be generated using the MATLAB function normrnd($\mu , \sigma ,\ldots$) which generates random variables with a normal distribution. The collection of random forces comprise $\varvec{\varGamma }_{B}$. The range of these forces are shown in Fig. 11a–c. These randomly fluctuating discontinuous functions slow numerical integration so each random variable is held constant during a single integration step. The random variable is updated at the beginning of each step and is held constant regardless of how many adaptive steps are required to reach the end of an integration time step.

Fig. 11

Range of random forces acting on atom: a $B_{A1_{1}}$, b $B_{A1_{2}}$, and c $\bar{r}_O \ B_{A1_{3}}$