ARTICLES

Sci. Rep. Kanazawa Univ.
Vol. 63, pp. 15–28, 2019

Molecular dynamics study of free energy profile fordissociation of ligand from CA I active site
MS A RWANSYAH∗1,2 , Koichi KODAMA1 , Isman K URNIAWAN1,3 , Tatsuki K ATAOKA1 ,Kazutomo K AWAGUCHI1 , and Hidemi NAGAO1
1 Division of Mathematical and Physical Sciences, Graduate School of Natural Science and Technology,

Kanazawa University, Kanazawa 920-1192, Japan
2 School of Chemistry, University of Cokroaminoto Palopo, Palopo, Indonesia
3 School of Computing, Telkom University, Jl. Telekomunikasi, Terusan Buah Batu, Bandung, Indonesia(Received November 21, 2018 and accepted in revised form March 10, 2019)

Abstract We investigate the binding/dissociation process of ligand molecule from carbonic anhydrase (CA) I carbonic anhydrase (CA) I enzyme by using all-atom molecular
dynamics simulation. The force field parameters of zinc ion in the CA I active site are estimated by quantum chemical calculations and are summarized in this paper. The free energy
profile for binding/dissociation process of ligand from CA I active site is calculated by thethermodynamic integration combined with the all-atom molecular dynamics simulation. The
binding free energy as a function of the distance between the center of mass positions of CA
I active site and the ligand molecule is estimated. The radial distribution function of the
CA I-ligand complex is calculated from the trajectory of all-atom molecular dynamics (MD)simulation. We estimate the free energy surface from the radial distribution function. We
can obtain the bond constant of the equilibrium state from the value of the free energy surface. We discuss the binding/dissociation process of ligand molecule by calculating the freeenergy profile to know the stability of the CA I-ligand complex with some thermodynamic
properties such as the binding free energy, the equilibrium state of the free energy surface
and so on.
Keywords. CA I active site, ligand molecule, force field parameters, free energy profile,
free energy surface.

1

Introduction

Carbonic anhydrases (CAs) family is a ubiquitous zinc enzyme which can be isolated from archaea, prokaryotes, and eukaryotes [1]. In higher vertebrates including human, 14 different CA
isozymes are found in several subcellular localizations such as CAI-III and CA VII found in cytosols, CA IV, CA IX, CA XII, and CA XIV located in membrane, CA V found in mitochondria,
and CA VI secreted in saliva [1, 2]. These enzymes involve the biological fundamental roles suchas catalyzes of the hydration of carbon dioxide to bicarbonate which is essential to regulate the pH
levels in cells, biosynthetic reaction and electrolyte secretion in several tissues [1–4].
∗ Corresponding author

Email: arwan@wriron1.s.kanazawa-u.ac.jp.

16

Molecular dynamics study of free energy profile for dissociation of ligand from CA I active site

The structures of carbonic anhydrases have been investigated by X-ray analysis [5–9]. All
isozymes of carbonic anhydrase contain a zinc ion in the enzyme active site which is essential forcatalytic activity. Three histidine residues coordinate to Zn ion through their imidazole nitrogenatoms. A water molecule occupies the fourth coordinated position in the active site at an acidic
pH (< 7) or a hydroxyl ion at higher pH with tetrahedral coordination geometry through protontransfer reaction from water binding to the zinc ion in the CAs active site [10, 11]. The selected
inhibitor can be applied to change the position of the water molecule or hydroxyl group to inhibit
the biological process in cells. The sulfonamide group of a heterocyclic structure can also be
applied in the primary application as the most potent inhibitors that bind to the zinc ion throughthe deprotonated nitrogen. The sulfonamide group is implemented to have an anticonvulsant,antiglaucoma, anticancer, and antiurolithic [4].
Almost every CAs family relates to some diseases such as glaucoma, diabetes, cancer, epilepsy,
and etc [12, 13]. The several ligand molecules to inhibit CAs activity are the critical target for thetherapeutics against many diseases. Therefore, understanding the interaction between ligand andCAs including Zn ion is crucial. The several computational studies on CAs have been recently performed [14–16]. The interaction of some antiepileptic drugs with CA active center has been quantum mechanically investigated to elucidate the relative orientations of CA II-inhibitor complex using a high-level calculation to understanding on the mechanism of inhibitor action [14]. All-atom
molecular dynamics (MD) simulation of CA enzyme has also been carried out for CA II enzyme
in complex with ligand (acetazolamide) to obtain a theoretical understanding of metalloproteininhibitor complex [15]. This paper shows the quantitative insight into the binding interaction ofprotein-ligand complex and represents the key role in Zn ion in the CAs active site. The theoreticalstudies on the standard free energy of protein-ligand binding have been investigated to understandthe free energy change of the ligand molecules in the binding/dissociation process [15–19].
In this paper, we investigate the binding free energy between CA I enzyme and the ligand
molecule. The force field parameters of the zinc ion in the CA I active site is estimated by quantum
chemical calculations. The free energy profile for binding/dissociation process of ligand fromCA I active site is calculated by using thermodynamic integration combined with all-atom MD
simulation according to the procedure provided in the previous work [20, 21]. We discuss thebinding/dissociation process of ligand molecule by calculating the free energy profile to know thestability of CA I-ligand complex related to some thermodynamic properties such as the binding
free energy, the equilibrium state of the free energy surface and so on.

2

Method and model

In this section, we present a model using in this simulation conditions for MD simulation. The
force field parameters in the CA I active site are estimated, and the free energy profile is introduced.
The free energy profile as a function of the distance between the center of mass positions of
CA I active site and the ligand molecule, rcm , is calculated by all-atom MD simulation with thethermodynamic integration method presented in the previous work [19, 22].

2.1

Model for molecular dynamics simulation

In this study, the structure of CA I in complex with the ligand molecule (4-carboxyethyl benzenesulfonamide ethyl ester) is obtained from the X-ray crystallographic structure in the protein databank solved by Srivastava, et al. (PDB ID: 2NN7) and is applied to initial structure for the simulation of ligand binding/dissociation process [23]. In the CA I enzyme, the total number of residues

MS A RWANSYAH, K. KODAMA, I. K URNIAWAN, T. K ATAOKA, K. K AWAGUCHI, H. NAGAO

17

(b)

(a)

Figure 1: (a) The structure of CA I in complex with the ligand, Zn2+ , and water molecules are placed intoa cubic box. CA I is represented by cartoon model. Ligand and Zn2+ ion are represented by VDW modeland silver spheres, respectively. (b) The cluster model of CA I active site in Zn2+ ion which is tetrahedrally
coordinated with three histidine residues (His 94, 96, 119) and ligand molecule.

is 260 which consists of a zinc ion, 10 helices and 18 strands of the beta sheet with the total number of atoms being 2029. The ligand is a derivative of benzenesulfonamide CA inhibitor whichcontains 17 atoms. The structure of CA I enzyme and ligand immersed by water molecules areshown in Fig. 1(a).

2.2Force field parameters in metal

The cluster model of CA I active site consists of a zinc ion and three histidine residues (His 94,
96, 119) which are tetrahedrally coordinated with a ligand molecule as shown in Fig. 1(b). The
CA I active site in complex with ligand is calculated by the quantum chemical methods to know
the interaction effect of the stability from the interaction of zinc ion with histidine residues andligand molecule.
The force field parameters of zinc ion in the CA I active site are not available in the ZincAMBER force field (ZAFF) database developed by Peters et al [24]. Therefore, we evaluate theforce field parameters of zinc ion in the CA I active site by calculating the potential energy surface(PES) as a function of bond distance and bond angle. The bond distance and angle potential areprovided as follows:
V (r, K r) = K r(r − rc)2 ,(2.1)
V (θ , K θ ) = K θ (θ − θ c)2 ,(2.2)

where K r and K θ are the force constants of distance and angle respectively, and rc and θ c are equilibrium distance and angle, respectively. The structure of CA I is optimized by freezing the heavy
atoms. The potential energy surface of bond distance and bond angle of CA I is analyzed by using
the B3LYP method with the 6-31G* basis set with the increment of 0.02 and 2◦ , respectively. The
range of increment is determined to avoid the inharmonic effect of distance and angle according to
the suggestion in Ref. [25, 26]. The optimized structure and potential energy surface calculations
for distance and angle are performed by Gaussian 09 packages [27]. Then, we convert the quantum
chemical calculation of the distance and angle of zinc-N(His 94, 96, 119) motif in the CA I active
site into the Amber force field parameters by using Metal Center Parameter Builder (MCPB) [28].
Our force field parameters have not been submitted to the AMBER parameter database.

18

2.3

Molecular dynamics study of free energy profile for dissociation of ligand from CA I active site

Simulation condition

All-atom molecular dynamics simulation is performed on CA I enzyme and ligand molecule. The
TIP3P water model [29] with 9808 water molecules inserted in a 70 × 69 × 77 Å3 periodic box
then, 3 Cl – ions are added to neutralize the whole system. The total number of atoms in the system
is 33,425. The AMBER14 force field is used for the protein molecule [30], and general AMBER
force field (GAFF) is applied to determine the force field parameters of ligand molecule. The
electrostatic interactions are computed by using the Particle Mesh Ewald (PME) algorithm. The
switching cutoff distance is 10 Å. The SHAKE algorithm is used to constrain the bond distance ofthe hydrogen atoms. A time step of 2 fs is applied in all simulations.
The energy minimization of the system is carried out with the constraint on the position of theheavy atoms of CA I, ligand and Cl – ion. Then, we perform the energy minimization without any
constraint. The system is simulated on NV T -constant simulation for 60 ps where the temperatureis gradually increased from 0 to 300 K. The temperature and pressure of the system are kept
at 300 K and 1 atm by using the Langevin thermostat and isotropic position scaling algorithm,
respectively. We equilibrate the system with the NPT ensemble for 50 ns without the harmonic
positional constraint. After the whole system reaches the equilibrium state, the production run
with a constraint on the distance between the centers of mass of CA I active site and the ligand
is performed. The SHAKE method is used to constraint the distance. We prepared 19 distanceswith the increment of 0.5 Å: r= 6-15 Å. Each simulation is performed for 5 ns with NPT -constantMD simulation to obtain the mean force < F(r0 ) > from the MD trajectory. All MD simulations
are carried out by Amber 16 packages [30]. The analysis of the trajectories of MD simulation is
performed by CPPTRAJ tools [31].

2.4

Root-mean-square deviation

We calculate the root-mean-square deviation (RMSD) of Cα atoms of CA I enzyme and the heavy
elements of ligand molecule according to the following equation:

RMSD(t 1) =

1

2
N

1 ∑
mi k ri t 1 − rre f ,i k2 ,
M i=1

(2.3)

where mi is the mass of atom i, N is the total number of Cα atoms of CA I enzyme and the heavy
elements of ligand, M is the total mass of Cα atoms in CA I and the heavy elements of ligand
molecule, ri (t) is the position of the Cα atom of protein and C, N, O, and S atoms of ligand at thetime t, and rre f ,i is the position of ith atom in the X-ray structure.

2.5

Free energy calculation

We calculate the ligand dissociation process from the CA I active site by generating the free energy
profile ∆G(r) as a function of the distance, rcm , between the centers of mass of ligand and CA I
active site as follows [32, 33]:
∆G(r) = G(r) − G(r0)dG(r0 ) 0
)dr
dr0
r0

Z r

=

(

Z r

=

<
r0

∂U
>r0 dr,0
∂ r0

(2.4)

MS A RWANSYAH, K. KODAMA, I. K URNIAWAN, T. K ATAOKA, K. K AWAGUCHI, H. NAGAO

19

where U is the potential energy of the whole system, r0 is a reference distance, G(r) is the freeenergy of the whole system, and < ... >0r represents the isothermal-isobaric ensemble average,
where the distance rcm is constrained to r0 . The free energy profile is evaluated by [33–35] asfollows :∆G(r) = −

Z r
r0

0

F(r ) =

< F(r0 ) >r0 dr0 ,


mpro
mlig
Flig −
Fpro · n,mpro + mlig
mpro + mlig

(2.5)

where F(r0 ) is the mean force at rcm , Flig and Fpro are forces at rcm , mpro and mlig are the total
masses of the CA I active site and ligand, respectively, and n is the unit vector from the centerof mass of the CA I active site to that of the ligand molecule. Then, ∆G(r) can be calculatedby [20, 36] as follows:
Nw
1∆G(r) ∼= − ∑ (< F(r0 ) >r0 =ri + < F(r0 ) >r0 =ri−1 )(ri − ri−1 ).
i=1 2

(2.6)

The free energy profile of the ligand dissociation as a function of the distance between the centersof mass of the CA I active site and the ligand molecule can be calculated using < F(r0 ) > fromtrajectories of MD simulations.

3

Results and Discussion

In order to investigate the dissociation process of the ligand molecule from the CA I active site,
we estimate the thermodynamic properties of ligand binding/dissociation which can be obtained
from the free energy profile as a function of rcm . In this section, we determine the force fieldparameters of the zinc ion in the CA I active site for MD simulation and then estimate the freeenergy profile of the ligand molecule in the binding/dissociation process by the thermodynamic
integration method. Finally, we discuss the dynamic of the complex in equilibrium state related to
the free energy surface estimated from the radial distribution function.

3.1

Force field parameters

The structure of model cluster as shown in Fig. 2 in CA I active is optimized to obtain a good
relative spatial arrangement of its structure. The residues (His 94, 96, 119) of CA I active site incomplex with ligand are connected to the Zn center represented by the tetrahedral coordination
structure. The schematic diagram of the charge distribution of the CA I active site is also shownin Fig. 3.
The potential energy surfaces (PES) of bond and of angle of CA I active site are shown in
Fig. 4. The equilibrium bond distances of Zn-N(His94), Zn-N(His96) and Zn-N(His119) arearound 1.999, 2.024 and 2.018 Å respectively as shown in Fig. 4(a). The fitted force constant(Kq) for Zn-N(His94), Zn-N(His96), and Zn-N(His119) is 93.419, 96.543, and 96.135 kcal mol−1
Å−2 , respectively. The PES of the equilibrium angle of X-Zn-Y (X,Y) = (N(His94), N(His96)),(N(His94), N(His119)), (N(His96), N(His119)) was 108.531◦ , 112.021◦ , 108.172◦ , respectively(Fig. 4(b)). The fitted force constant (Kq) for (His94)N-Zn-N(His96), (His96)N-Zn-N(His119),and (His94)N-Zn-N(His119) is 72.351, 77.055, and 71.863 kcal mol−1 rad−2 , respectively. Table 1 summarizes the bond and angle parameters of zinc ion in the CA I active site from the20

Molecular dynamics study of free energy profile for dissociation of ligand from CA I active site

Figure 2: The optimized structure of CA I active site in complex with ligand molecule in which the zinc
ion is located at the center of a deep, silver spheres model, tetrahedrally coordinated with three histidine
residues (His 94, 96, 119) and ligand molecule.

Figure 3: Charge distribution around the CA I active site.

experimental values and quantum chemical calculations. The experimental values of bond andangle parameters are estimated from the mean values (with standard deviations) of the eight X-ray
crystallographic structures of CA I sharing the Zinc-N(His 94, 96, 119) motif. From Table 1, we
find that the values of the equilibrium distances for Zn-N(His94), Zn-N(His96) and Zn-N(His119)
calculated by quantum chemical method was similar to those obtained by the experimental results,
and the angle of (His94)N-Zn-N(His96) is the most favorable parameter with similar values of theexperimental results. The bond and angle force field parameters of zinc ion in the CA I active siteare evaluated by quantum chemical calculation according to the calculation procedure in Ref. [37]and [38].

3.2

Molecular dynamics simulation

The root-mean-square deviation (RMSD) value defined in Equation (2.3) is evaluated to check thestability of the CA I in complex with the ligand molecule from the initial modeling structure by
X-ray analysis. Fig. 5(a) shows the time evolution of the RMSD value which is obtained from
the MD simulation of the CA I in complex with the ligand, and the simulated system almost
becomes the equilibrium state at 2 ns. From Fig. 5(b), we find that the average distance of rcm

MS A RWANSYAH, K. KODAMA, I. K URNIAWAN, T. K ATAOKA, K. K AWAGUCHI, H. NAGAO

21

(b)

(a)

Figure 4: The potential energy surface (PES) (a) PES for bond distance between Zn and X (X=N(His94,
His96, His94)), (b) PES for bond angle X-Zn-Y ((X,Y)= (N(His94), N(His96)), (N(His94), N(His119)),(N(His96), N(His119))

Table 1: Bond and angle parameters of zinc ion in the CA I active site. The standard deviations are shownin parentheses.

Bond DistanceZn-N(His94)
Zn-N(His96)
Zn-N(His119)
Bond Angle(His94)N-Zn-N(His96)(His94)N-Zn-N(His119)(His96)N-Zn-N(His119)
Zn-N(His94)-Cβ
Zn-N(His94)-Cγ
Zn-N(His96)-Cβ
Zn-N(His96)-Cγ
Zn-N(His119)-Cβ
Zn-N(His119)-Cγ

rexperimenta(Å)
2.040 (0.051)
2.125 (0.083)
2.063 (0.054)
θ experimenta(◦ )107.175 (2.539)
110.958 (4.859)
102.174 (5.296)
121.166 (4.450)
130.514 (5.240)
138.397 (6.745)
112.348 (6.839)
131.993 (2.190)
117.769 (1.975)

rQC
(Å)
1.999
2.024
2.018
θ QC
(◦ )
108.531112.021108.172
126.848126.571128.203125.474
127.997124.874

Kq(kcal mol−1 Å−2 )93.41996.54396.135Kθ
(kcal mol−1 rad−2 )72.35177.05571.863198.938198.838204.398
96.542201.740
191.955

a The experimental data obtained from the X-ray crystallographic structures of CA I sharing the Zinc-

N(His 94, 96, 119) motif in the protein data bank (PDB ID: 2NN7 [23], 2NN1 [23], 2NMX [23],1AZM [39], 1BZM [39], 3LXE [40], 4WUQ [41], 6EVR [42].)for the last 3.0 ns is 8.58 Å with the standard deviation of 0.29 Å. Fig. 6 shows the fluctuationof the distances and angles in the CA I active site. We obtain that the distances of zinc-N(His94, 96, 119) motif are 1.95± 0.05, 2.00± 0.05, and 2.01± 0.06 Å and the angles of X-Zn-Y((X,Y)= (N(His94), N(His96)), (N(His94), N(His119)), (N(His96), N(His119)) are 107.73± 3.14,
113.64± 2.76, 103.40± 2.94, respectively. These results are consistent with those of the distances
and angles from quantum chemical calculation in Table 1. We assume that the distances and angles
of zinc-N(His 94, 96, 119) motif in the CA I active site become stable during the 50 ns simulation.

22

Molecular dynamics study of free energy profile for dissociation of ligand from CA I active site
(b)

(a)

Figure 5: (a) The RMSD value of CA I-Ligand complex, (b) the distance between the center of mass ofligand and that of CA I active site during the simulation.
(a)

(b)
3

180

Zn−N(His94)
Zn−N(His96)
Zn−N(His119)

2.8

160
Angle (Degree)

2.6
Distance (Å)(His94)N−Zn−N(His96)(His94)N−Zn−N(His119)(His96)N−Zn−N(His119)

140

2.4
2.2

120

2

100

1.8
80

1.6
0

10

20

30

Time (ns)

40

50

0

10

20

30

40

50

Time (ns)Figure 6: (a) The distance and (b) angle fluctuation in the CA I active site during the simulation.

3.3

Free energy profile

The mean force < F(r0 ) > derived in Equation (2.5) is estimated at each rcm values shown in
Fig. 7(a). The error bars represent the standard deviation which is obtained from the independent measurements of < F(r0 ) > for each 5-ns trajectory of MD simulation. The repulsive force
between CA I active site and ligand is represented as the positive value of < F(r0 ) >, while thenegative value represents an attractive force between CA I active site and ligand. The attractive
force between CA I active site and ligand molecule is observed at the distance of 9 Å ≤ rcm ≤ 15
Å, while the repulsive force at the distance of 6 Å rcm 8.5 Å. We assume that the distance of 15 Å isthe reference state, because the mean force at the position becomes small at around 5.4 × 10−11 N.
It implies that the interaction between CA I and ligand molecule almost vanishes at the reference
state. We observe that the error bar at rcm 14.5 Å is smaller than the other distances, and it suggeststhat the conformational change of ligand molecule has small fluctuation around the side chains of
CA I active site and, on the other hand, that the error bar at rcm 11 Å has large value because of
the large fluctuation of the ligand conformation in the side chains of CA I active site.
Fig. 7(b) shows the free energy profile as a function of rcm derived in Equation (2.6). Here,
we discuss the binding free energy from the reference state (rcm = 15 Å) and not from the standard
state of the binding free energy. The binding free energy from the reference state is calculated byusing thermodynamic integration method to elucidate the stabilization of the ligand molecule in thebinding pocket of CA I enzyme. The binding free energy as shown in Ref. [43] has been calculated

23

MS A RWANSYAH, K. KODAMA, I. K URNIAWAN, T. K ATAOKA, K. K AWAGUCHI, H. NAGAO
(b)

(a)

20

5
4

0

2

∆G(r) (kcal/mol)<F(r)> (10−9 N)

3
1
0
−1
−2
−3

−20
−40
−60

−4
−5
5

6

7

8

9

10

11

12

13

14

15

16

−80

5

Distance (Å)

6

7

8

9

10

11

12

13

14

15

Distance (Å)

Figure 7: (a) The mean force < F(r0 ) > and (b) free energy profile G(r) as a function of the distance
between the center of mass of ligand and that of CA I active site. The error bars represent the standard
deviation for all the trajectories of MD simulation

by a method to calculate the binding free energy with reference to the standard state. The standard
chemical potential of each molecule is defined as a hypothetical standard state where the moleculein the system occupies at standard concentration in the reaction solvent. From Fig. 5(b), we canapproximately estimate the binding free energy of the complex from the reference state of -53.11±
14.17 kcal/mol and the free energy reaches the minimum at rcm 8.5 Å. This result is consistent withthat of the equilibrium MD, in which the average value is rcm 8.58 Å. Meanwhile, the calculated
binding free energy (-53.11 kcal/mol) is underestimated compared to Ref. [15]. However, forthe case of the binding free energy of protein/ligand complex, the energy could be around -53.11
kcal/mol or more [19, 44, 45]. Unfortunately, we cannot validate at present the estimated binding
free energy because there is no corresponding experimental result available for the binding freeenergy along the distance between CA I/ ligand complex. However, we believe that our results
contribute to giving an insight for estimation of the free energy profile as the function of thebinding distance.

3.4

Free energy surface around the binding state

In the previous subsection, we have estimated the free energy profile corresponding to the binding/dissociation process between the CA I active site and the ligand as shown in Fig. 7(b) and haveestimated the binding free energy between the molecules. However, the standard deviation of freeenergy is not so small at the equilibrium state. Therefore, it is not easy to know the detail of behavior of the dynamics of CA I and ligand complex at the equilibrium state of the association from
Fig. 7(b). For the purpose, we try to investigate the dynamics of the complex from the viewpointof another physical approach.
Fig. 8 shows the radial distribution function (RDF) as a function of rcm . The radial distributionfunction (RDF) g(rcm ) can be calculated from the trajectory of all-atom MD simulation for thelast 3 ns. The value of the first peak is represented as g(rcm ) = ρ(rcm )/ρ max, where ρ(rcm ) andρ max represent the number of density at the distance rcm and the highest peak value, respectively.
We find that the position of the peak of the radial distribution function is 8.58 Å. The value of rcm
in the X-ray analysis is 8.58 Å. This result is the similar value obtained from RDF calculation asshown in Fig. 8. The free energy surface can be estimated as G(rcm ) = kBT lng(rcm ), where kB, T ,and g(rcm ) correspond to the Boltzmann constant, temperature and the radial distribution function
for the distance between the centers of mass of CA I active site and ligand molecule, respectively.

24

Molecular dynamics study of free energy profile for dissociation of ligand from CA I active site
1.5

1

1.2
∆G(r) (kcal/mol)

1.2

g(r)

0.8
0.6
0.4

0.6
0.3

0.2
0

0.9

4

5

6

7

8

9

10

11

12

13

14

0

8

8.2

8.4

8.6

8.8

9

9.2

Distance (Å)

Distance (Å)

Figure 8: The radial distribution function (RDF)

Figure 9: The free energy profile obtained fromof the distance between the center of mass of ligand and that of CA I active site.

RDF data in Figure 6.

Fig. 9 shows the free energy surface obtained from the radial distribution function, where thevalue of the free energy is calculated from the equilibrium state, and the horizontal axis represents
the distance between the center of mass of CA I active site and that of ligand molecule. The
lowest free energy value along the distance is observed at the distance of 8.58 Å and we can findthat the minimum distance is similar to that obtained by thermodynamic integration calculations.
From Fig. 9, we can find the harmonicity of the free energy around the equilibrium state and canestimate the bond constant from the equilibrium state of the free energy surface as 271.72 kcal
mol−1 Å−2 at the distance of 8.58 Å. The bond potential is represented as F(rcm ) = K(rcm − r0)2 ,where K and r0 are the bond constant and equilibrium distance, respectively.

3.5

Binding Free Energy based on MM/PBSA calculation

The binding free energy between CA I and ligand from the reference state has been calculated
in the previous subsection. We have also calculated the radial distribution function as a functionof rcm . The free energy surface around the equilibrium state has been determined from the RDF
value. However, we can not obtain the binding free energy from the standard state. Therefore,
we investigate the binding free energy of the protein-ligand complex in water solvent by estimating the gas-phase interaction energy EMM , solvation free energy ESolv , and entropy −T ∆S of thecomplex. These energies represent as ∆Gbind = ∆EMM + ∆GSolvPBSA − T ∆S. The gas-phase interaction energy between CA I and ligand is given by ∆EMM = ∆Eele + ∆Evdw , where ∆Eele and∆Evdw are the van der Waals and the electrostatic energies, respectively. The solvation free energy ∆GPBSA
Solv is estimated by using continuum approach (Poisson-Boltzmann/surface area), i.e.,
GSolvPBSA = GSolvPB +GSA . The contribution energies as shown in the equation above are calculated
Solv
by using MM/PBSA program and Normal mode analysis (Nmode) [46] in the Amber 16 packages. All the binding free energies of CA I/ligand complexes are calculated by using a single
configuration obtained from the procedure described in the section above (simulation condition).
Table 2 lists the contribution of the molecular mechanics and solvation energies of the CA I
enzyme obtained from the several PDB files with some different ligands. From Table 2, we obtainthat the binding free energy for all CA I/ligand complexes are good aggrement with the experimental results. For the fourth of the binding free energies, the van der Waals energies are notsignificantly different, which implies they all have good hydrophobic interactions. However, wefind that the electrostatic energies (∆Eele ) are different with values of -119.25, -166.54, -129.35,-130.79 kcal/mol. We conclude that the electrostatic term influences for determining the different
binding orientations.

25

MS A RWANSYAH, K. KODAMA, I. K URNIAWAN, T. K ATAOKA, K. K AWAGUCHI, H. NAGAO

Table 2: The binding free energy and the contribution of each energy term. The standard deviations areshown in parentheses. Units are in kcal/mol.
PDB ID
2NN7
2NN1
2NMX
1AZM

Total Charge
-3.00
-2.00
-1.00
-1.00

∆Evdw
-16.93-10.60-14.66-11.51∆Eele-119.25-166.54-129.35-130.79

∆GPB
Solv
89.25
130.5697.70
103.57

∆GSA
Solv
-3.21
-2.64
-2.93
-2.49

EMM
-136.18-177.14-144.01-143.87∆GPBSA
Solv
86.04
127.9294.77
101.07−T ∆S-41.50-41.90-41.72-34.72

Gcalc
bind
-8.64
-7.32
-7.52
-8.08

Exp

Gbind
-9.57 a-8.43 a-7.61 a-8.29 bThe experimental values for the binding free energy of CA I/ligand complexes are obtained from the Ref. [23]a and Ref. [39]b .

4

Summary

In this study, we have performed the all-atom MD simulation of the CA I complex with the ligand
molecule. We have presented a simple cluster model derived from the structure by X-ray analysis
to obtain the force field of the zinc ion in the CA I active site. The force field parameters related
to the zinc ion with MD simulation has been summarized. The free energy profile of the binding/dissociation process of ligand from the CA I enzyme has been estimated by some integration
methods. We have discussed the free energy surface of the CA I enzyme with the ligand in relation
to the radial distribution function of the distance between the centers of mass of CA I enzyme andthe ligand molecule.
In the estimation of free energy profile, the mean force acting on between the CA enzyme andthe ligand has been estimated at each distance between the center of the molecule from the MD
simulation. We have found that the attractive force can be observed at the longer distance thanthe equilibrium distance and that the repulsive force can be observed in the shorter distance thanthe equilibrium distance. From the mean force, we have calculated the free energy profile. The
results of free energy profile suggest that the equilibrium distance derived from MD simulations isa good agreement with the experimental distance, and we have estimated the binding free energyof -53.11± 14.17 kcal/mol. Meanwhile, we calculate the free energy surface between the CA I
enzyme and the ligand by the radial distribution function estimated from the MD trajectories. We
have found the harmonicity of the free energy surface around the equilibrium distance between
two molecules with the bond constant of 271.72 kcal mol−1 Å−2 at the equilibrium point at the
distance of 8.58 Å.
We have also calculated the binding free energy based on MM/PBSA calculation to estimate thecontribution of the molecular mechanics and solvation energies of several CA I/ligand complexes.
We observe that the electrostatic energy is the predominant impact for determining the binding
free energy of the complexes.

Acknowledgements
The first author expresses thanks to Ministry of Research, Technology and Higher education(KEMRISTEKDIKTI) and Indonesia Endowment Fund for Education (LPDP) through a BUDILN Scholarship Scheme for financial support during a Ph.D study at Kanazawa University, Japan.

References[1] C. T. Supuran, F. Briganti, S. Tilli, W. Chegwidden, A. Scozzafava, “Carbonic anhydrase inhibitors: sulfonamides
as antitumor agents?” Bioorg. Med. Chem 9 (3) (2001) 703–714.
[2] C. T. Supuran, A. Scozzafava, A. Casini, “Carbonic Anhydrase Inhibitors,” Med. Res. Rev 23 (2) (2003) 146–189.

26

Molecular dynamics study of free energy profile for dissociation of ligand from CA I active site

[3] C. T. Supuran, A. Scozzafava, “Applications of carbonic anhydrase inhibitors and activators in therapy,” Expert
Opin. Ther. Pat. 12 (2) (2002) 217–242.
[4] C. T. Supuran, A. Scozzafava, “Carbonic anhydrase inhibitors and their therapeutic potential,” Expert Opin. Ther.
Pat. 10 (5) (2000) 575–600.
[5] S. K. Nair, J. F. Krebs, D. W. Christianson, C. A. Fierke, “Structural basis of inhibitor affinity to variants of humancarbonic anhydrase II,” Biochemistry 34 (12) (1995) 3981–3989.
[6] D. A. Whittington, A. Waheed, B. Ulmasov, G. N. Shah, J. H. Grubb, W. S. Sly, D. W. Christianson, “Crystal structure of the dimeric extracellular domain of human carbonic anhydrase xii, a bitopic membrane protein
overexpressed in certain cancer tumor cells”, Proc. Natl. Acad. Sci. U. S. A. 98 (17) (2001) 9545–9550.
[7] M. Ferraroni, S. Tilli, F. Briganti, W. R. Chegwidden, C. T. Supuran, K. E. Wiebauer, R. E. Tashian, A. Scozzafava,“Crystal structure of a zinc-activated variant of human carbonic anhydrase I, CA I Michigan 1: Evidence for asecond zinc binding site involving arginine coordination,” Biochemistry 41 (20) (2002) 6237–6244.
[8] J. Leitans, A. Kazaks, A. Balode, J. Ivanova, R. Zalubovskis, C. T. Supuran, K. Tars, “Efficient Expression and
Crystallization System of Cancer-Associated Carbonic Anhydrase Isoform IX,” J. Med. Chem. 58 (22) (2015)
9004–9009.
[9] P. Řezáčová, D. Borek, S. F. Moy, A. Joachimiak, Z. Otwinowski, “Crystal structure and putative function ofsmall Toprim domain-containing protein from Bacillus stearothermophilus,” Proteins: Structure, Function, andBioinformatics 70 (2) (2008) 311–319. doi:10.1002/prot.21511.
[10] W. S. Sly, P. Y. Hu, “Human carbonic anhydrases and carbonic anhydrase deficiencies,” Annu. Rev. Biochem. 64(1995) 375–401.
[11] S. Pastorekova, S. Parkkila, J. Pastorek, C. T. Supuran, “Carbonic anhydrases: Current state of the art, therapeutic
applications and future prospects,” J. Enzyme Inhib. Med. Chem. 19 (3) (2004) 199–229.
[12] A. Scozzafava, L. Menabuoni, F. Mincione, F. Briganti, G. Mincione, C. T. Supuran, “Carbonic anhydrase inhibitors. synthesis of water-soluble, topically effective, intraocular pressure-lowering aromatic/heterocyclic sulfonamides containing cationic or anionic moieties: is the tail more important than the ring?,” J. Med. Chem.
42 (14) (1999) 2641–2650.
[13] C. T. Supuran, A. Scozzafava, “Carbonic anhydrases as targets for medicinal chemistry,” Bioorg. Med. Chem.
15 (13) (2007) 4336–4350.
[14] M. Ghiasi, S. Kamalinahad, M. Arabieh, M. Zahedi, “Carbonic anhydrase inhibitors: A quantum mechanicalstudy of interaction between some antiepileptic drugs with active center of carbonic anhydrase enzyme,” Comput.
Theor. Chem. 992 (2012) 59–69.
[15] T. O. Wambo, L. Y. Chen, S. F. McHardy, A. T. Tsin, “Molecular dynamics study of human carbonic anhydrase
II in complex with Zn2+and acetazolamide on the basis of all-atom force field simulations,” Biophys. Chem.
214-215 (2016) 54–60.
[16] Y. Deng, B. Roux, “Computations of Standard Binding Free Energies with Molecular Dynamics Simulations,” J.
Phys. Chem. B 113 (8) (2009) 2234–2246. doi:10.1021/jp807701h.
[17] S. Doudou, N. A. Burton, R. H. Henchman, “Standard free energy of binding from a one-dimensional potentialof mean force,” J. Chem. Theory Comput. 5 (4) (2009) 909–918.
[18] I. J. General, “A Note on the Standard State’s Binding Free Energy,” J. Chem. Theory Comput. 6 (8) (2010)
2520–2524.
[19] K. Kawaguchi, H. Saito, H. Nagao, S. Okazaki, “Molecular dynamics study on the free energy profile for dissociation of adp from n-terminal domain of hsp90,” Chem. Phys. Lett. 588 (2013) 226–230.
[20] K. Kawaguchi, H. Saito, H. Nagao, “Decomposition analysis of free energy profile for Hsp90-ADP association,”
Mol. Simul. 42 (11) (2016) 896–901.
[21] S. Nakagawa, I. Kurniawan, K. Kodama, M. S. Arwansyah, K. Kawaguchi, H. Nagao, “Theoretical study oninteraction of cytochrome f and plastocyanin complex by a simple coarse-grained model with molecular crowding
effect,” Mol. Phys. 116 (5-6) (2018) 666–677.
[22] K. Kawaguchi, S. Nakagawa, S. Kinoshita, M. Wada, H. Saito, H. Nagao, “A simple coarse-grained model forinteracting protein complex,” Molecular Physics 115 (5) (2017) 587–597.
[23] D. K. Srivastava, K. M. Jude, A. L. Banerjee, M. Haldar, S. Manokaran, J. Kooren, S. Mallik, D. W. Christianson,“Structural analysis of charge discrimination in the binding of inhibitors to human carbonic anhydrases i and ii,”
J. Am. Chem. Soc. 129 (17) (2007) 5528–5537.

MS A RWANSYAH, K. KODAMA, I. K URNIAWAN, T. K ATAOKA, K. K AWAGUCHI, H. NAGAO

27[24] M. B. Peters, Y. Yang, B. Wang, L. Füsti-Molnár, M. N. Weaver, K. M. Merz, “Structural Survey of ZincContaining Proteins and Development of the Zinc AMBER Force Field (ZAFF),” Journal of Chemical Theory
and Computation 6 (9) (2010) 2935–2947. doi:10.1021/ct1002626.
[25] C. Liu, B. Zhang, Y. Zhu, M. Tang, “Evaluations of AMBER force field parameters by MINA approach forcopper-based nucleases,” Struct. Chem. 27 (5) (2016) 1449–1464.
[26] Y. Zhu, Y. Su, X. Li, Y. Wang, G. Chen, “Evaluation of Amber force field parameters for copper(II) withpyridylmethyl-amine and benzimidazolylmethyl-amine ligands: A quantum chemical study,” Chem. Phys. Lett.
455 (4-6) (2008) 354–360.
[27] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone,
B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino,
G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima,Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J.
Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C.
Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken,
C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W.
Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich,
A. D. Daniels, O. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski, D. J. Fox, “Gaussian 09, Revision E.01.”[28] P. Li, K. M. Merz, “Mcpb.py: A python based metal center parameter builder,” Journal of Chemical Informationand Modeling 56 (4) (2016) 599–604.
[29] W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, M. L. Klein, “Comparison of simple potentialfunctions for simulating liquid water,” J. Chem. Phys. 79 (2) (1983) 926–935.
[30] R. Salomon-Ferrer, D. A. Case, R. C. Walker, “An overview of the amber biomolecular simulation package,”
WIREs. Comput. Mol. Sci. 3 (2) (2013) 198–210.
[31] D. R. Roe, T. E. Cheatham, “Ptraj and cpptraj: Software for processing and analysis of molecular dynamics
trajectory data,” J. Chem. Theory Comput. 9 (7) (2013) 3084–3095.
[32] M. Zacharias, T. P. Straatsma, J. A. McCammon, F. A. Quiocho, “Inversion of Receptor Binding Preferencesby Mutagenesis: Free Energy Thermodynamic Integration Studies on Sugar Binding to l-Arabinose Binding
Proteins,” Biochemistry 32 (29) (1993) 7428–7434.
[33] Y. Andoh, K. Oono, S. Okazaki, I. Hatta, “A molecular dynamics study of the lateral free energy profile of a pairof cholesterol molecules as a function of their distance in phospholipid bilayers,” J. Chem. Phys. (2012) 155104.
[34] K. Fujimoto, N. Yoshii, S. Okazaki, “Free energy profiles for penetration of methane and water molecules into
spherical sodium dodecyl sulfate micelles obtained using the thermodynamic integration method combined withmolecular dynamics calculations,” J. Chem. Phys. (2012) 094902.
[35] S. Kawada, K. Fujimoto, N. Yoshii, S. Okazaki, “Molecular dynamics study of the potential of mean force of SDS
aggregates,” J. Chem. Phys. (2017) 084903.
[36] K. Kawaguchi, S. Nakagawa, I. Kurniawan, K. Kodama, M. S. Arwansyah, H. Nagao, “A coarse-grained modelof the effective interaction for charged amino acid residues and its application to formation of gcn4-pli tetramer,”Mol. Phys. 116 (5-6) (2018) 649–657.
[37] T. Shuku, K. Sugimori, A. Sugiyama, H. Nagao, T. Sakurai, K. Nishikawa, “Molecular orbital analysis of active
site of oxidized azurin: Dependency of electronic properties on molecular structure,” Polyhedron 24 (16-17)(2005) 2665–2670.
[38] I. Kurniawan, T. Matsui, K. Kawaguchi, H. Nagao, “Theoretical studies on association/dissociation process ofplastocyanin and cytochrome f in photosynthesis,” J. Phys.: Conf. Ser. (in press).
[39] A. Casini, J. Y. Winum, J. L. Montero, A. Scozzafava, C. T. Supuran, “Carbonic anhydrase inhibitors: Inhibitionof cytosolic isozymes I and II with sulfamide derivatives,” Bioorganic and Medicinal Chemistry Letters 13 (5)(2003) 837–840.
[40] V. Alterio, S. M. Monti, E. Truppo, C. Pedone, C. T. Supuran, G. De Simone, “The first example of a significant
active site conformational rearrangement in a carbonic anhydrase-inhibitor adduct: the carbonic anhydrase itopiramate complex,” Org. Biomol. Chem. 8 (2010) 3528–3533.
[41] A. Zubrienė, J. Smirnovienė, A. Smirnov, V. Morkũnaitė, V. Michailovienė, J. Jachno, V. Juozapaitienė,P. Norvaišas, E. Manakova, S. Gražulis, D. Matulis, “Intrinsic thermodynamics of 4-substituted-2,3,5,6tetrafluorobenzenesulfonamide binding to carbonic anhydrases by isothermal titration calorimetry,” Biophysical
Chemistry 205 (2015) 51–65. doi:https://doi.org/10.1016/j.bpc.2015.05.009.

28

Molecular dynamics study of free energy profile for dissociation of ligand from CA I active site[42] N. Chiaramonte, M. Romanelli, E. Teodori, C. Supuran, “Amino Acids as Building Blocks for Carbonic Anhydrase Inhibitors,” Metabolites 8 (2) (2018) 36.
[43] M. K. Gilson, J. A. Given, B. L. Bush, J. A. McCammon, “The statistical-thermodynamic basis for computationof binding affinities: A critical review,” Biophysical Journal 72 (3) (1997) 1047–1069.
[44] S. K. Choubey, D. Prabhu, M. Nachiappan, J. Biswal, J. Jeyakanthan, “Molecular modeling, dynamics studies
and density functional theory approaches to identify potential inhibitors of sirt4 protein from homo sapiens : anovel target for the treatment of type 2 diabetes,” Journal of Biomolecular Structure and Dynamics 35 (15) (2017)
3316–3329.
[45] K. A. Brown, Y. Zou, D. Shirvanyants, J. Zhang, S. Samanta, P. K. Mantravadi, N. V. Dokholyan, A. Deiters,“Light-cleavable rapamycin dimer as an optical trigger for protein dimerization,” Chem. Commun. 51 (2015)
5702–5705.
[46] B. R. Miller, T. D. Mcgee, J. M. Swails, N. Homeyer, H. Gohlke, A. E. Roitberg, “MMPBSA . py : An Efficient
Program for End-State Free Energy Calculations,” J. Chem. Theory Comput. 8 (9) (2012) 3314–3321.