Detachment of HCO3– from the Active Site of Carbonic Anhydrase

Aug 13, 2018 - Detachment of HCO3– from the Active Site of Carbonic Anhydrase: ... results, suggesting the great potential of ANN for re-engineering...
0 downloads 0 Views 924KB Size
Subscriber access provided by Kaohsiung Medical University

C: Physical Processes in Nanomaterials and Nanostructures 3–

Detachment of HCO from the Active Site of Carbonic Anhydrase: Molecular Dynamics Simulation and Machine Learning Gong Chen, Diannan Lu, Jianzhong Wu, and Zheng Liu J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.8b05298 • Publication Date (Web): 13 Aug 2018 Downloaded from http://pubs.acs.org on August 18, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Detachment of HCO3– from the Active Site of Carbonic Anhydrase: Molecular Dynamics Simulation and Machine Learning Gong Chen1, Diannan Lu1*, Jianzhong Wu2, Zheng Liu1*1 1Key

Lab of Industrial Biocatalysis, Ministry of Education, China, Dept. of Chemical Engineering, Tsinghua University, Beijing 100084, China,

2

1

Dept. of Chemical and Environmental Engineering, University of California, Riverside, CA 92521, USA

Email: [email protected]; [email protected] 1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract While tremendous efforts have been devoted, for many years both experimentally and theoretically, to exploring the potential of carbonic anhydrase (CA) as an industrial catalyst for CO2 sequestration, little is known about the kinetics of HCO3– detachment from the active site of the enzyme, an essential step affecting the performance of the CO2 hydration process. Using molecular dynamics simulation and the Markov-state model, we established a global kinetic profile for HCO3– detachment from CA. Using human carbonic anhydrase II (hCA-II) as the model enzyme, we calculated the commitors of all Markov states for HCO3– detachment through the analysis of transition paths, and found that HCO3– prefers to diffuse toward the bulk solution once it dissociates with the binding sites. We also established the kinetic map for HCO3– detachment using molecular simulation data with the Markov-state model and identified six distinct pathways for HCO3– diffusion towards the bulk solution. These pathways were determined according to the electrostatic interaction gradient and the solvation gradient of HCO3– produced by amino acid residues of hCA-II and the surrounding water molecules, respectively. To explore the relationship of HCO3– residence time at different binding sites as a function of steric hindrance, local charge density and hydrophobicity, we proposed an artificial neural network (ANN) algorithm to correlate HCO3– occupancy probability, which reflects binding strength in different binding sites, with the molecular properties of hCA-II. The predicted HCO3– occupancy probability agrees perfectly with the simulation results, suggesting the great potential of ANN for re-engineering carbonic anhydrase as well as other enzymes toward tailored applications.

2

ACS Paragon Plus Environment

Page 2 of 32

Page 3 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Introduction The diffusion of substrates and/or products during enzymatic catalysis has profound effects on the reaction outcome and kinetics1, 2. Carbonic anhydrase (CA) catalyzes transformation of CO2 to HCO3– in water at an extremely high rate and thus is of great interest for CO2 sequestration3-13. Although a good understanding of CO2 diffusion form the bulk phase to the active site of CA has been established by both experimental observation3, 14-16 and molecular simulation17-34, a comprehensive view of HCO3– diffusion from the active site of the enzyme to the bulk phase, however, is yet to be established. The numerical analysis of the relationship between the kinetics of HCO3– detachment and the catalysis performance of CA shows that a smaller rate constant for HCO3– detachment would lead to a drastic reduction of the hydration outcome. This indicates that HCO3– detachment has significant effects on the kinetics of CA-catalyzed CO2 hydration (see Supporting Information for details). It has been shown that3 HCO3– dissociation from CA follows two steps, i.e., the breakage of Zn2+-HCO3– coordination bond and HCO3– detachment from the active site. For the breakage of coordination bonds, Merz et al.35 studied various binding modes between hCA-II and HCO3– and showed that Zn2+-HCO3– binding can be described in terms of the Lindskog and Lipscomb modes. The exchange between these two modes results in destabilization of the Zn2+-HCO3– coordination bond by Thr199 and, subsequently, HCO3– cleavage and its replacement by a water molecule. Regarding HCO3– detachment from the active site of CA, Piazzeta et al36 examined the initial detachment of HCO3– and H2NCOHN- in the active site of hCA-II by molecular dynamics (MD) simulation. The potential mean force (PMF) for HCO3– shows that HCO3– is more readily to detach from the active site than H2NCOHN–. Despite all these efforts, a global kinetic profile of HCO3– diffusion pathway is not yet available, which is of essential importance to understand and, eventually, re-engineer CA for targeted 3

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

applications. Molecular dynamics (MD) simulations have been previously used to investigate binding, diffusion and detachment of small molecules anchored in protein binding pockets37-40. Simulation yields useful information such as the free-energy changes due to the binding processes. However, large fluctuations in protein configurations and the lack of statistical significance in short-time MD sampling make it difficult to calculate the kinetic constants for the entire diffusion process (exceeding 100 nanoseconds). The Markov-state model (MSM)41 are suitable for studying molecular processes with high energy barriers. Multiple Markov states are introduced to eschew direct sampling of energy barriers thus enable calculation of the transfer rates between metastable states and the evolution of the time-evolved probability in each Markov states. Previously, we used MSM to investigate the equilibrated properties and the kinetics of complex events, e.g. the diffusion of CO2 towards to the active site of CA42, 43. Artificial neural network (ANN) has been widely used in chemical and biochemical engineering fields as a powerful tool to explore and predict complex relationship between the target and featured variables with excellent accuracy and universality. Both experimental and simulated data can be used for machine learning44, 45. A combination of computational chemistry and machine learning algorithms has also been applied to screening of catalysts46, 47, and design and evaluation of gas adsorbents48 and drug molecules49. With its proven effectiveness in extracting potential features of complex events through different hidden layers, we come to the idea of using the nonlinear multi-layered complex network model, in conjunction with MD and MSM, to probe the complex interaction of HCO3– with CA. The correlation intrinsically underpins the release kinetics and pathway of HCO3– association with CA. 4

ACS Paragon Plus Environment

Page 4 of 32

Page 5 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

In the present study, we use MD and MSM to analyze the detachment pathway of HCO3– dissociation from the active site of hCA-II, from which we obtain the diffusion kinetics network, the driving force and the diffusion rate constant. In addition, we employ ANN to correlate the equilibrium occupancy of each Markov state with the local environment properties in terms of the probability of cavity opening, the local charge density, and the local hydrophobicity index. ANN allows us to reproduce the occupancy probability obtained from the Markov-state model. With these correlations we are able to establish a global kinetic model for HCO3– detachment with molecular details.

Models and Methods Molecular Model The molecular structure of hCA-II was obtained from the PDB database (PDB code: 2CBA50). We used H++51 to determine the protonation state for each side chain of amino acids in hCA-II at pH=9 because CA is utilized as an additive to promote the rate of CO2 absorption/desorption in the amine-based post combustion CO2 capture process (at pH 9~108, 9). The configurations of zinc-histidine and zinc-water complex were the same as those in our previous work43. The parameters of HCO3– ions were generated from CGenff52 developed for the CHARMM force field, and TIP3P53 was chosen as the water model, other molecules were all described with CHARMM36 force field54. In all simulations, we used the genbox command in the Gromacs software to randomly insert 100 HCO3–, 99 Na+ and 15406 water molecules around the protein centered in the simulation box. The simulation box size was set as 8 nm×8 nm×8 nm. Simulation Methods All MD simulations were performed using the Gromacs 4.5.4 software55. V-rescale thermostat and Parrinello-Rahman barostat were used to keep the system at 293 K and 1 atm, with the coupling time of 5

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

0.2 ps and 2.0 ps, respectively. A switch method was adopted to smooth the Lennard-Jones potential at the cut-off distance of 1.2 nm. The electrostatic interaction was calculated through the Particle Mesh Ewald (PME) method with the grid size of 0.16 nm. The integration time step in MD was 2 fs, simulation data were collected every 1000 steps. Periodic boundary conditions were applied to all 3 dimensions. In each MD run, we first minimized the system energy to make sure that the maximum force in the system was less than 1500 kJ mol-1 nm-1. Then the pre-equilibrium simulation was conducted for 10 ps with the constraints on the backbone carbon atoms of hCA-II. Another pre-equilibrium simulation was followed for 10 ps without any constraints. Finally, MD simulation was performed for 600 ns for data collection. To assure that the sampling of the 600 ns trajectory is sufficient for profiling alternative pathways, we have performed 5 additional MD simulations with different initial positions of 100 HCO3– with the duration of 50 ns. The radical distribution function of CA was compared with that obtained from the sampling of the 600 ns trajectory. Figure S3 shows a perfect agreement, indicating that sampling of the 600 ns trajectory is sufficient for generating alternative pathways. (See supporting information for details). All snapshots were generated by VMD 1.9.156. Markov State Model In our previous work, we used a Markov state model (MSM) to study CO2 diffusion to the active site of hCA-II 43. The trajectory was split into 100 ones such that each short trajectory spans 600 ns with only 1 HCO3–. Before performing the clustering algorithm, we expressed all trajectories to position hCA-II at the center of the simulation box. Moreover, at each moment hCA-II was fitted to its initial position by minimizing the RMSD between the instantaneous conformation and the initial conformation. In the meantime, the HCO3– molecule underwent the same translation and rotation 6

ACS Paragon Plus Environment

Page 6 of 32

Page 7 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

moves with the instantaneous configuration of hCA-II. We used 2 Å as the cut-off radius for the cluster analysis, and only considered the HCO3– molecule within 6 Å from the α-C atom of hCA-II to avoid getting some Markov states too far away from CA. We identified 101 non-overlapping Markov states with a radius of 2 Å and one bulk phase Markov state through the cluster analysis process (by g_cluster command in Gromacs55). To identify the hCA-II binding site for HCO3–, we utilized the crystal structure of hCA-II, PDB code 2VVB57, that contains the binding of HCO3– in the active-site pocket. The HCO3– binding site in 2CBA could be identified by best fitting of the 2VVB structure to 2CBA (Figure 1). We set the initial position of HCO3– diffusion determined by the 2VVB structure as Markov state 1, which is the starting point of HCO3– detachment, while the other Markov states, which are obtained by analyzing the MD trajectories with the clustering algorithm, are ranked according to the distance from Markov state 1. Figure 1 shows the good fitting of 2VVB to 2CBA, indicating that the position of HCO3– in 2VVB could be utilized in our MSM. After the division of the MD trajectories, we obtained a total of 102 Markov states in hCA-II and one bulk-phase Markov state. The amino-acid composition of each Markov state is defined as all amino acids whose distance from the Markov state centroid is less than or equal to 2.0 Å plus their van der Waals radii.

7

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 32

Figure 1. Fitting results for the Markov states. To get the precise position of the HCO3– binding site in the active site of CA, the 2CBA structure was best fitted to the 2VVB structure. The binding sites, which represent the Markov states of 2CBA, are shown in transparent green spheres. The structures of 2CBA(cyan) and 2VVB(blue) were both drawn in NewCartoon mode, while HCO3– in the 2VVB system was shown in CPK mode. We used the method described previously43 to determine the probability of cavity opening for each Markov state. An opening status refers to the condition that the Markov state can maintain a cavity with radius of 2.0 Å that does not overlap with the van der Waals radii of atoms of all surrounding amino acids. Once the composition of the surrounding Markov states was determined, we defined the state hydrophobicity, that is, the sum of the hydrophobicity indices of all amino acids involved in the Markov state. Each amino-acid hydrophobicity index is judged according to Hessa's method58. Because HCO3– is negatively charged, the occupancy of HCO3– may have a high correlation with the local charges of amino acids around the Markov state itself. Therefore, the overall charge of all the amino acids within 5 Å of the centroid of the Markov state is assumed as the local charge density of the Markov state. The master equation is adopted to evolve the time-dependent probabilities in all Markov states59:

  d  (t )  K   (t ) dt

(1)

 where K is a rate transfer matrix,  (t ) is a vector composed of the occupation probabilities of all Markov states. Each element in the transfer rate matrix, kij, represents the rate of transition from the Markov state j to the Markov state i; it satisfies the following two conditions:

kij  0 i  j

k

ij

(2)

0

(3)

i

The transfer rate from Markov state j to i is calculated from the MD trajectory with the following 8

ACS Paragon Plus Environment

Page 9 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

equations59:

N ij

kij 

, j  i

N

Tj



l 1,l  j

N lj

k jj  

1 Tj

(4)

(5)

where Nij represents the number of transitions from Markov state j to i, and Tj is the ensemble average residence time of HCO3– molecules at Markov state j. We assume that the probabilities of the bulk-phase Markov state can be described by the first-order kinetic equation pbulk (t ) 

k1 [1  exp( ( k1 +k 2 [HCO 3 ])t )]  k1 +k 2 [HCO 3 ]

(6)

where k1 is the first-order rate constant for HCO3– diffusion into the bulk phase, and k2 is the second-order rate constant for HCO3– diffusion into the active state of hCA-II, [HCO3–] is the HCO3– concentration set as 1/(NAꞏVbox) in our MSM analysis, NA is the Avogadro constant, and Vbox is the volume of the simulation box. Transition Path Theory

In order to identify the main path of HCO3– diffusion from the active state, we used the transition path theory60 based on the Markov transition rate matrix. For HCO3– diffusion, we divided the Markov states into three categories: 1) the diffusion initiation state (state B): HCO3– Markov state 1; 2) intermediate diffusion states (state I): all other Markov states (2-102) in hCA-II except for Markov state 1; 3) and the diffusion endpoint (state S): the bulk-phase Markov state. According to the Markov state model and the transition path theory, all commitors of state I could be obtained by solving a series 9

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 32

of linear equations60, 61:

 iI

B

(i)kij   kij  0, j  I

(7)

iB

The reactive flux between any two Markov states can be calculated from J j i  kij p eqj [B (i )  B ( j )] B (i )  B ( j ) .

(8)

Artificial Neural Network Algorithm

We used Artificial Neural Network (ANN) to correlate HCO3– occupancy to the local properties of hCA-II represented by the cavity opening probability, hydrophobicity, and local charge density of the Markov state. The ANN algorithm in our study (as shown in Fig. S2), also called multi-layer perceptron, is a nonlinear network optimization model62. We added a nonlinear activation function to the signal transfer: Relu function, which enables the model to learn nonlinear data. The learning rate of ANN was set to 1E-4. Adam function63, which was a stochastic gradient-based optimizer proposed by Kingma et al., was adopted as the solver function; alpha factor (L2 regular penalty function parameter) was set to 0.2, the number of hidden layer was 3 with 50 nodes (neurons) per layer, and a maximum number of iterations for the solver was set to 20,000. The Relu function is given by:

Relu(x)  max(x, 0)

(9)

Before training, we preprocessed the input labeled data by using the StandardScaler method in Scikit-learn62, which made the average value of each labeled data equals to 0 and the variance of each labeled data is normalized to 1.0. We used the Python-based machine learning library scikit-learn62 to build the ANN algorithm. 90% of the data from molecular simulation and from the Markov-state model were used as the training set for our ANN model. The remaining 10% of the data were used as 10

ACS Paragon Plus Environment

Page 11 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

the testing set.

Results and Discussion Validation of the Markov state model

Through MD simulation and the MSM analysis, we can depict the diffusion paths of HCO3– from the CA binding site to the bulk phase. To establish the Markov-state model, first we applied different lag times, i.e., 2 ps, 4 ps, 10 ps, 20 ps, and 40 ps to build the transition-rate matrix K. Evolved probability of each Markov state was then calculated according to the master equation (1) and the initial probability vector. It was found that the time-dependent probability evolution curves of different Markov states have different patterns. As shown in Fig. 2, the evolution curve tends to be consistent for a typical Markov state in different lag times, indicating that HCO3– diffusion is Markovian, and the lag time converges in the interval of 2–40 ps. Therefore, the lag time of 10 ps was chosen in the following analysis.

11

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 2. Time-evolved probabilities for Markov states 1, 2, 22 and the bulk phase state under different lag times (2 ps, 4 ps, 10 ps and 40 ps) Figure 2 shows three different types of evolution curves. The evolution curve for Markov state 1 is presented in Fig. 2(a); it is referred to as the “quick reduction type” as the time-evolved probability of Markov state 1 monotonically decreases and approaches 0 within 1 ns, i.e., HCO3– ions could be quickly removed from Markov state 1. As shown in Fig. 2(b), the time-evolved probability for Markov state 2 belongs to "pathway type", indicating that the time-evolved probability of HCO3– is not monotonous. Markov state 2 might be a binding site in the diffusion pathways of HCO3– from the active site to the bulky solution. Fig.2(c) shows that the evolution curve of Markov state 22 belongs to the “growth type”. The time-evolved probability for Markov state 22 monotonically increases, which means the HCO3– occupancy in this state is continuously accumulating. This suggests that Markov 12

ACS Paragon Plus Environment

Page 12 of 32

Page 13 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

state 22 has a certain affinity towards HCO3–. As shown in Fig. 2(d), the time-evolved probability for the bulk phase state also belongs to the “growth type”. Detachment track of HCO3-

(a) “Active site”

0

“Bulk phase state”

“intermediate state”

1

(c)

(b) 3

9 6

1 2 1.00 0.94 0.88 0.82 0.76 0.70

Figure 3. (a) Schematic of the transition path theory used in this work, (b) Commitors of Markov states in all hCA-II, (c) Mapping of the Markov states and their commitors in the hCA-II. The Markov states in the figure are represented by vdw balls (red spheres). The color of each vdw ball is correlated with the value of its commitor. In order to examine the detachment pathway of HCO3– from the active site of hCA-II, we calculated the commitors of all Markov states based on equation (7). Fig. 3(a) gives a schematic diagram for application of the transition path theory to the kinetics of HCO3– detachment from the active site, in which the commitor of the "active site" in the diffusion initiation state, Markov state 1, is set to 0, and the commitor of the destination, i.e., the bulk phase Markov state, is set to 1.0. The commitors for all other states could be calculated from Eq. (7). The Markov state with a larger commitor is closer to the end of the path in the context of “diffusion”. Fig. 3(b) shows the commitors for different Markov states. It indicates that the commitors of Markov states 2-102 ranged from 0.79 to 13

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 32

1.0, suggesting that these Markov states are remote from the active site and close to the bulk phase. Once the HCO3– ions move beyond Markov state 1, there is a substantial probability that they may diffuse into the bulk phase. On the other hand, the Markov states with small sequence numbers 2, 3 have smaller commitors than any other Markov states (except for Markov state 1), suggesting that these Markov states are closer to Markov state 1 than other intermediate Markov states in terms of “diffusion”. Fig. 3(c) showed the positions and commitors of Markov states in the hCA-II. Except for Markov states 1, 2, and 3, the commitors of other Markov states are close to 1.0. As shown in Fig. 3(c), Markov states 2 and 3, which have lower commitor values than any other Markov states except Markov state 1, are located in the natural hydrophobic pocket of hCA-II.

(a)

(b) 3

9

6

2

1 1

2 3 3

9 6

6

1 2

9

Figure 4 (a) Markov states 1,2,3,6 and 9 in the diffusion pathway of HCO3–, red arrows denote the reactive fluxes calculated from eqn (8), their directions represent the direction of reactive fluxes and their sizes are set according to the flux values, (b) time-evolved probability curves of Markov states 1, 2, 3, 6, 9 and bulk phase states on the diffusion path. After obtaining the commitors of all Markov states, we calculated the reactive fluxes between different Markov states according to equation (8) and identified the diffusion path of HCO3– in terms of these Markov states. As illustrated in Figure 4(a), the reactive fluxes of Markov states 2, 3, 6 and 9 to 14

ACS Paragon Plus Environment

Page 15 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

the bulk phase state accounted for 99.1% of the reactive flux from the Markov state 1 to the bulk phase state. Therefore, these Markov states form a diffusion pathway for HCO3– detachment; they are drawn in the hCA-II structure as shown in the upper panel of Figure 4(a). Figure 4(a) also shows the direction and the magnitude of the reactive fluxes in hCA-II, in which the magnitude is indicated by the width of arrowed cylinders (Figure 4(a), the bottom panel). These reactive fluxes clearly depict the diffusion pathway of HCO3–. Fig.4 (b) presents the time-evolved probability curves of Markov states 1, 2, 3, 6, 9 and the bulk phase state. It shows that the probability of Markov state 1 decreases continuously, and the probabilities of Markov states 2, 3, 6 and 9 exhibit bell-shape curves along the diffusion path, whereas the probability of the bulk phase state increases monotonically. The continuous change of HCO3– occupancy probability for each Markov state indicates that the transition path theory is fully consistent with the time-evolved occupancy probabilities obtained from the Markov-state model. To investigate HCO3– diffusion, we analyzed the kinetic network in the diffusion process and the corresponding snapshots of HCO3– in each Markov state. In Figure 5, we assign different areas for the Markov states according to their equilibrium occupancy probabilities. Here the bulk-phase state has the highest occupancy probability, whereas Markov state 1 has the lowest probability. Markov states 2,3,6,9 have greater occupancy probabilities than Markov state 1, indicating that their significant role in determining the HCO3– diffusion pathways. With this method, we are able to identify 4 distinct diffusion pathways: 1) 1→2→bulk phase state, 2) 1→3→bulk phase, 3) 1→3→ 6→bulk phase, and 4) 1→3→9→bulk phase. These multiple diffusion pathways indicate that HCO3– in Markov state 1, which is generated from the Zn2+ coordination structure, does not diffuse directly to the bulk phase. On the contrary, it needs to go through at least one intermediate state, i.e., Markov state 2, 3, 6 or 9.

15

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 5. Diagrammatic representation of the HCO3– diffusion network. Here each dot represents the Markov state in the diffusion pathway, the dot area corresponds to the relative size of the equilibrium occupancy probability in this Markov state, and the arrow represents the direction of transition with its size proportional to the reactive flux calculated according to equation (8). The composition of the amino acid residues in the Markov state, HCO3– ions present in this Markov state, and water molecules in the surrounding 5 Å are represented by a stick model, hydrogen bonds between water molecules and HCO3– ion are indicated by red dotted lines. Figure 5 presents both the direction and the magnitude of the reactive fluxes. The reactive flux from Markov state 1 to 2 is 1.3E-6, which is 4.8-fold higher than that from Markov state 1 to 3. Moreover, the reactive flux from Markov state 2 to 3, 1.1E-6, is also significantly higher than those 16

ACS Paragon Plus Environment

Page 16 of 32

Page 17 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

corresponding to other pathways, indicating the high probability of HCO3– diffusion from Markov state 2 to 3. and then the subsequent steps has more options, e.g., through Markov states 6 and 9, or directly from 3 to bulk solution. The direct diffusion to the bulk phase accounts for 53.9% of the reactive fluxes diffused from Markov state 3, while 21.6% and 24.5% are towards Markov states 6 and 9, respectively. Figure 5 also gives a snapshot for each Markov state occupied by HCO3–. These snapshots provide molecular insight on the local compositions of the Markov states. For example, for HCO3– in Markov 1, the hydroxyl group of HCO3– interacts with threonine (Thr263), and HCO3– forms 3 pairs of hydrogen bonds with water molecules in the surrounding within 5 Å. For HCO3– in Markov state 2, the surrounding amino acids include two histidines (His64, His94) and one asparagine (Asn62). HCO3– formed 4 pairs of hydrogen bonds with the surrounding water molecules. For HCO3– in Markov state 3, the surrounding amino acids include tryptophan (Trp5), histidine (His64), threonine (Thr200), and asparagine (Asn62), all of which are hydrophilic except tryptophan. HCO3– forms 6 pairs of hydrogen bonds with the surrounding water molecules. For HCO3– in Markov state 6, the surrounding amino acids consist of leucine (Leu60) and arginine (Arg58). This Markov state is less sterically hindered; HCO3– forms two pairs of hydrogen bonds with Arg58, and 8 pairs of hydrogen bonds with water molecules. For HCO3– in Markov state 9, it interacts with lysine (Lys170). In this case, the steric hindrance is insignificant; and HCO3– forms 8 pairs of hydrogen bonds with the surrounding water molecules.

17

ACS Paragon Plus Environment

The Journal of Physical Chemistry

-400

BC-(Pro+Water)

11 10

-420

9 -440

8

-460

7 6

-480 -500

# of Hydrogen Bonds

Nonboning Energy(KJ/mol)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 32

5 0

2

4 6 Cluster ID

8

4 10

Figure 6. The non-bonding energies (electrostatic + vdw) of HCO3– due to its interaction with the surrounding amino-acid residues and water molecules in Markov states 1, 2, 3, 6 and 9 along different diffusion pathways, and the number of hydrogen bonds between HCO3–and the surrounding water molecules. To analyze the driving force for HCO3– detachment from the active site, we calculated the non-bonding energies between HCO3– and amino acids from CA and water molecule when HCO3– is in Markov states 1, 2, 3, 6 and 9, respectively. Fig. 6 shows the number of hydrogen bonds between HCO3– and the surrounding water molecules. For HCO3– in Markov state 1, the non-bonding energy for its interaction with the surrounding amino acid residues and water molecules is –422.8 kJ/mol, which is the weakest point in the HCO3– diffusion path. The non-bonding energies of HCO3– at Markov states 2, 3 are –449.1 kJ/mol and –469.9 kJ/mol, respectively. The electrostatic interaction between HCO3– and the surrounding amino acid residues and water molecules dominates the total non-bonding energy. The number of hydrogen bonds between HCO3– and the surrounding water molecules is 6.4, 8.2 and 8.4 at Markov states 1, 2 and 3, respectively. Looking back to the aforementioned six diffusion pathways, 1) 1→2→3→bulk phase, 2) 1→2→ 3→6→bulk phase, 3) 1→2→3→9→bulk phase, 4) 1→3→bulk phase, 5) 1→3→6→bulk phase, 6)1 18

ACS Paragon Plus Environment

Page 19 of 32

→3→9→bulk phase, we find that HCO3– experiences a greater non-bonding energy with hCA-II

residues and forms more pairs of hydrogen bonds with the local water molecules as it diffuses to the bulk solution. Therefore, the solvation effects for HCO3– become stronger along the diffusion paths from the active site to the bulk solution. In other words, the electrostatic interaction gradients produced by the surrounding amino acids and solvation effects gradients induced by water molecules in different binding sites of the pathway serve as the major driving force for HCO3– detachment from the active site of hCA-II.

Occupied Probability (Bulk phase state)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

1.0 0.8 0.6 0.4 0.2

Fitting Occupied Probability

0.0 0.01

0.1

1

10

100

1000 10000

Time (ns) Figure 7. The evolution of the Markov state for HCO3– in the bulk phase. The solid line represents the fitting of the time-dependent probability curve with equation (6) To obtain the rate constant for HCO3– diffusion from the binding site to the bulk solution, we fit the time-evolved probability curve with equation (6), and the results are shown in Fig. 7. The proposed first-order rate equation could well describe the time-evolved probability. In the following, equation (10) represents the kinetics of HCO3– detachment from the active site and its diffusion to 19

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 32

the bulk phase in the CO2 hydration cycle in hCA-II. Here [hCA-II:HCO3– ] represents the physical binding state of HCO3– at the active site, and k1 and k2 are a pair of reversible diffusion rate constants, Kd is the physical dissociation constant. [hCA  II: HCO3 ]

k1 k2

[hCA  II]  [HCO3 ]

Kd 

(10)

k1 k2

(11)

Table 1 Kinetics and thermodynamic parameters of HCO3- diffusion from hCA-II active site

Parameters

values

k1 (s-1)

2.63×108

k2 (mM-1 s-1)

8.94×106

Kd (mM)

29.4

Table 1 presents the diffusion rate constants and the equilibrium constant in equations (10) and (11) obtained through the fitting. Compared with the rate constant for CO2 diffusion in hCA-II43, we found that HCO3– diffuses slower towards the bulk solution than CO2. Besides, HCO3– is more difficult to physically dissociate from the active site of hCA-II due to the stringer electrostatic interaction of HCO3– with hCA-II.

20

ACS Paragon Plus Environment

Page 21 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Correlation of HCO3- occupied probability with the local environment

Figure 8. Correlation between the HCO3- equilibrium occupancy of all Markov states with (a) the cavity open probability of the states, (b) the local charge density, (c) the hydrophobicity index, (d) the HCO3- average residence time of the states. According to the results shown in Figures 3-7, there are multiple pathways, with different driving forces and rate constants, for HCO3– diffusion from the active site of hCA-II to the bulk solution. ANN allows us to predict the dependence of the binding strength of HCO3– on hCA-II structure. Such knowledge is essential for the molecular reengineering of CA-II. First, we analyzed the relations between equilibrium occupancy of HCO3– with the cavity open probability of the Markov state, the local charge density, the hydrophobicity index, and the average residence time of HCO3– in the Markov states. As shown in Fig. 8 (a), the equilibrium occupancy is relatively low for Markov states with a cavity open probability less than 20%, all of them are below 21

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

0.002. As for those Markov states with the cavity open probability above 20%, their occupancy probabilities varied by a large margin, indicating that weak steric hindrance is a necessary but insufficient condition for high HCO3- equilibrium occupancy in the Markov state. Fig. 8 (b) correlates the equilibrium occupancy in the Markov state with the local charge density, i.e., the overall charges on all atoms of the protein within 5 Å from the geometry center of the Markov state. Fig. 8 (b) shows there is a weak but positive relation of occupancy of the Markov state with the local charge density. In other words, HCO3– prefers to occupy Markov states with a positive charge density. Fig. 8 (c) shows that the equilibrium occupancy of HCO3– in the Markov state exhibits a weak but positive correlation with the hydrophobicity of the Markov state, i.e., HCO3– tends to occupy a Markov state with a large hydrophobicity index, i.e., a more hydrophilic Markov state. It is noted that, as shown by Fig. 8 (d), a higher HCO3– equilibrium occupancy is obtained from MSM in response to a longer HCO3– average residence obtained from MD trajectory. However, there also exist some Markov states with relatively low equilibrium occupancy but a large number of average residence frames in the MD data.

22

ACS Paragon Plus Environment

Page 22 of 32

Page 23 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Prediction of HCO3- occupied probability by linear regression and ANN

Figure 9. Comparison between the simulation data and the model predictions of HCO3– occupancy on each Markov state: (a) training and learning results of a linear fitting model based on the number of average residence frames, (b) training and testing results of the ANN model based on the hCA-II cavity open probability, the hydrophobicity index, and the local charge density Table 2. Average unsigned error and root mean square error in the linear regression model model

Data sets

Average unsigned error

Root mean square error

Linear

Training sets

6.16×10-4

8.23×10-4

Regression

Testing sets

4.73×10-4

7.29×10-4

Training sets

7.64×10-6

7.33×10-5

Testing sets

4.14×10-4

5.11×10-4

ANN

The binding sites and their relative binding strengths (occupancy probabilities) play an important role in determining the catalytic mechanism of hCA-II for the conversions of small molecules such as HCO3–. Fig. 8 indicates that the HCO3– residence time, obtained from the MD data, is positively correlated with the HCO3– equilibrium occupancy. Therefore, we proposed a linear regression for the prediction of the HCO3– occupancy in each Markov state. 23

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

According to the data set of the HCO3– equilibrium occupancy in terms of its corresponding HCO3– residence time shown in Fig. 8 (d), 90% of the data were used as training set for a linear regression model, the remaining 10% of the data were used as the testing set. Figure 9 (a) and Table 2 compare the HCO3– occupancy probabilities of both training set and testing set obtained from the Markov state model and from the linear regression. The results indicate that the prediction in both training and test sets were unsatisfactory and the linear regression model based on the HCO3– average residence frames either underestimates or overestimates the HCO3– equilibrium occupancy in each Markov state, suggesting that the average residence time calculated from MD is insufficient to predict the equilibrium occupancy probability. Artificial neural network (ANN) is useful to describe non-linear events and could extract different features in the input data through the addition of hidden layers. In this work, we use ANN to explore the complex dynamic interaction of HCO3– with hCA-II and with the surrounding water molecules. The features for ANN training were selected according to our understanding of the given event. For HCO3–detachment from the active site of CA, the detachment path and rate are mainly determined by HCO3–interaction with the surrounding atoms. These interactions include steric hindrance, and electrostatic and hydrophobic forces. The steric effect is a major factor affecting the entrance and exit of HCO3–, and the charge and hydrophobic forces determine the moving rate of the HCO3– molecule. Therefore, the criteria to justify the ANN features are: a) the fundamental characteristics of the protein, b) the relevance between a feature with the occupancy probability, c) data accessibility. Although the relevance between the residence time and occupancy is stronger than other three features we choose, it is not chosen as an ANN feature because it not a fundamental character of CA and a long-time MD trajectory is needed to obtain the data. Instead, we selected three features of 24

ACS Paragon Plus Environment

Page 24 of 32

Page 25 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

hCA-II as the input: the cavity open probability, hydrophobicity, index and local charge density. These features allow us to predict the occupancy probabilities of HCO3– in different HCO3– binding sites. As shown in Figure 8 for the linear regression of those features with the occupancy probability, all those selected ANN features have certain relevance with the occupancy probability (binding affinity). Our ANN model used three hidden layers with 50 neurons per layer with RELU, i.e., equation (9) as the activation function. The Adam function63 was selected as the gradient solver. The maximum number of iterations was 20,000, the learning rate was 1E-4, and the convergence accuracy was 1E-4. The input data were normalized and preprocessed with the StandardScaler method before upload to the neural network. We use 90% of the data as the training sets and the remaining 10% as the testing sets. As shown in Fig. 9 (b) and Table 2, the prediction of the ANN sets for training sets is satisfactory, i.e., the predicted HCO3– occupancy is almost identical to that obtained from MSM. The even distribution of the red points (testing sets) near the diagonal line suggests that the proposed ANN can extract key features in the training data through the hidden layer and the activation function. It is noted that the ANN may lead to over fitting as indicated by the different AUE and RMS values of the training set and the testing set. To prevent over fitting, we will try more advanced ANN models in the future work. While ANN offers a powerful tool for a rapid prediction of the binding strength, it yields no insights into the mechanism of enzyme kinetics because all features and the target are accommodated in the hidden layer as a black box. Nevertheless, ANN would allow us to predict the binding strength of a substrate and/or the kinetics for the release of a product as a function of the molecular structure of an enzymes under specific conditions such as different pHs, salt 25

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

concentrations, and temperatures.

Conclusion We propose a computational protocol that combines molecular dynamics simulation, the Markov state model, and artificial neutron network (ANN) to predict the kinetics of HCO3– detachment from the active site of hCA-II. We show that evolution of the Markov states for HCO3– diffusion to the bulk can be categorized into three types: "rapid reduction type", "pathway type" and "growth type". Using the transition path theory, we identify six distinctive HCO3– detachment pathways connecting the active site of hCA-II with the bulk solution. We find that the driving force for HCO3– detachment is affiliated with its electrostatic interactions with the surrounding amino acids and with the solvation effects induced by the local water molecules in different metastable states. ANN is able to correlate the binding strengths of HCO3– at different binding sites identified by the Markov diffusion network in terms of the steric hindrance, the local charge density, the local hydrophobicity index, and the HCO3– resident time. The ANN prediction of HCO3– occupancy probability reproduces perfectly that from the molecular simulation. The above results illustrate the strength of the computing protocol in establishing a global kinetic profile, with molecular insight on enzymatic reaction steps, and moreover, highlight the potential in its application to molecular reengineering of enzymes for customized applications.

Author Information Corresponding Author * E-mail: [email protected]; [email protected] Notes

The authors declare no competing financial interest. 26

ACS Paragon Plus Environment

Page 26 of 32

Page 27 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Supporting Information Numerical calculation on the catalytic performance of CA; Molecular model of HCO3–; Validation of sampling efficiency; Equilibrium occupancy probability of each Markov state; Schematic of ANN in this work

Acknowledgements This work is supported by the Chinese National Natural Science Foundation under grant No. 21520102008. The numerical calculations were performed at the “Exploration 100” platform supported by Tsinghua National Laboratory for Information Science and Technology.

References 1. Watt, E. D.; Shimada, H.; Kovrigin, E. L.; Loria, J. P., The mechanism of rate-limiting motions in enzyme function. Proceedings of the National Academy of Sciences 2007, 104 (29), 11981-11986. 2. Chen, G.; Kong, X.; Zhu, J.; Lu, D.; Liu, Z., How ABA block polymers activate cytochrome c in toluene: molecular dynamics simulation and experimental observation. Phys. Chem. Chem. Phys. 2015, 17 (16), 10708-10714. 3. Krishnamurthy, V. M.; Kaufman, G. K.; Urbach, A. R.; Gitlin, I.; Gudiksen, K. L.; Weibel, D. B.; Whitesides, G. M., Carbonic anhydrase as a model for biophysical and physical-organic studies of proteins and protein-ligand binding. Chemical reviews 2008, 108 (3), 946-1051. 4. Ji, X.; Su, Z.; Wang, P.; Ma, G.; Zhang, S., Tethering of Nicotinamide Adenine Dinucleotide Inside Hollow Nanofibers for High-Yield Synthesis of Methanol from Carbon Dioxide Catalyzed by Coencapsulated Multienzymes. ACS Nano 2015, 9 (4), 4600-4610. 5. Chen, X.; Wang, Y.; Wang, P., Peptide Induced Affinity Binding of Carbonic Anhydrase to Carbon Nanotubes. Langmuir 2015. 6. Merle, G.; Fradette, S.; Madore, E.; Barralet, J. E., Electropolymerized Carbonic Anhydrase Immobilization for Carbon Dioxide Capture. Langmuir 2014, 30 (23), 6915-6919. 7. Jo, B. H.; Seo, J. H.; Yang, Y. J.; Baek, K.; Choi, Y. S.; Pack, S. P.; Oh, S. H.; Cha, H. J., Bioinspired Silica Nanocomposite with Autoencapsulated Carbonic Anhydrase as a Robust Biocatalyst for CO2 Sequestration. ACS Catalysis 2014, 4 (12), 4332-4340. 8. Alvizo, O.; Nguyen, L. J.; Savile, C. K.; Bresson, J. A.; Lakhapatri, S. L.; Solis, E. O.; Fox, R. J.; Broering, J. M.; Benoit, M. R.; Zimmerman, S. A., Directed evolution of an ultrastable carbonic anhydrase for highly efficient carbon capture from flue gas. Proceedings of the National Academy of Sciences 2014, 111 (46), 16436-16441. 9. Vinoba, M.; Bhagiyalakshmi, M.; Grace, A. N.; Kim, D. H.; Yoon, Y.; Nam, S. C.; Baek, I. H.; Jeong, S. K., Carbonic Anhydrase Promotes the Absorption Rate of CO2 in Post-Combustion Processes. The Journal of Physical Chemistry B 2013, 117 (18), 5683-5690. 27

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

10. Sahoo, P. C.; Jang, Y.-N.; Lee, S.-W., Immobilization of carbonic anhydrase and an artificial Zn (II) complex on a magnetic support for biomimetic carbon dioxide sequestration. Journal of Molecular Catalysis B: Enzymatic 2012, 82, 37-45. 11. Zhang, S.; Zhang, Z.; Lu, Y.; Rostam-Abadi, M.; Jones, A., Activity and stability of immobilized carbonic anhydrase for promoting CO2 absorption into a carbonate solution for post-combustion CO2 capture. Bioresource Technology 2011, 102 (22), 10194-10201. 12. Vinoba, M.; Lim, K. S.; Lee, S. H.; Jeong, S. K.; Alagar, M., Immobilization of human carbonic anhydrase on gold nanoparticles assembled onto amine/thiol-functionalized mesoporous SBA-15 for biomimetic sequestration of CO2. Langmuir 2011, 27 (10), 6227-6234. 13. Vinoba, M.; Bhagiyalakshmi, M.; Jeong, S. K.; Yoon, Y., II; Nam, S. C., Capture and Sequestration of CO2 by Human Carbonic Anhydrase Covalently Immobilized onto Amine-Functionalized SBA-15. The Journal of Physical Chemistry C 2011, 115 (41), 20209-20216. 14. Bertini, I.; Luchinat, C.; Monnanni, R.; Roelens, S.; Moratal, J., Interaction of carbon dioxide and copper (II) carbonic anhydrase. Journal of the American Chemical Society 1987, 109 (25), 7855-7856. 15. Eriksson, A. E.; Kylsten, P. M.; Jones, T. A.; Liljas, A., Crystallographic studies of inhibitor binding sites in human carbonic anhydrase II: a pentacoordinated binding of the SCN− ion to the zinc at high pH. Proteins: Structure, Function, and Bioinformatics 1988, 4 (4), 283-293. 16. Fierke, C. A.; Calderone, T. L.; Krebs, J. F., Functional consequences of engineering the hydrophobic pocket of carbonic anhydrase II. Biochemistry 1991, 30 (46), 11054-11063. 17. Mikulski, R.; West, D.; Sippel, K. H.; Avvaru, B. S.; Aggarwal, M.; Tu, C.; McKenna, R.; Silverman, D. N., Water Networks in Fast Proton Transfer during Catalysis by Human Carbonic Anhydrase II. Biochemistry 2013, 52 (1), 125-131. 18. Maupin, C. M.; Castillo, N.; Taraphder, S.; Tu, C.; McKenna, R.; Silverman, D. N.; Voth, G. A., Chemical Rescue of Enzymes: Proton Transfer in Mutants of Human Carbonic Anhydrase II. Journal of the American Chemical Society 2011, 133 (16), 6223-6234. 19. Kaila, V. R. I.; Hummer, G., Energetics and dynamics of proton transfer reactions along short water wires. Phys. Chem. Chem. Phys. 2011, 13 (29), 13207-13215. 20. Hakkim, V.; Subramanian, V., Role of Second Coordination Sphere Amino Acid Residues on the Proton Transfer Mechanism of Human Carbonic Anhydrase II (HCA II). The Journal of Physical Chemistry A 2010, 114 (30), 7952-7959. 21. Roy, A.; Taraphder, S., Transition Path Sampling Study of the Conformational Fluctuation of His-64 in Human Carbonic Anhydrase II. The Journal of Physical Chemistry B 2009, 113 (37), 12555-12564. 22. Maupin, C. M.; McKenna, R.; Silverman, D. N.; Voth, G. A., Elucidation of the Proton Transport Mechanism in Human Carbonic Anhydrase II. Journal of the American Chemical Society 2009, 131 (22), 7598-7608. 23. Krishnamurthy, V. M.; Kaufman, G. K.; Urbach, A. R.; Gitlin, I.; Gudiksen, K. L.; Weibel, D. B.; Whitesides, G. M., Carbonic Anhydrase as a Model for Biophysical and Physical-Organic Studies of Proteins and Protein−Ligand Binding. Chemical Reviews 2008, 108 (3), 946-1051. 24. Swanson, J. M. J.; Maupin, C. M.; Chen, H.; Petersen, M. K.; Xu, J.; Wu, Y.; Voth, G. A., Proton Solvation and Transport in Aqueous and Biomolecular Systems:  Insights from Computer Simulations. The Journal of Physical Chemistry B 2007, 111 (17), 4300-4314. 25. Riccardi, D.; Cui, Q., pKa Analysis for the Zinc-Bound Water in Human Carbonic Anhydrase II:  28

ACS Paragon Plus Environment

Page 28 of 32

Page 29 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Benchmark for “Multiscale” QM/MM Simulations and Mechanistic Implications†. The Journal of Physical Chemistry A 2007, 111 (26), 5703-5711. 26. Riccardi, D.; König, P.; Prat-Resina, X.; Yu, H.; Elstner, M.; Frauenheim, T.; Cui, Q., “Proton Holes” in Long-Range Proton Transfer Reactions in Solution and Enzymes:  A Theoretical Analysis. Journal of the American Chemical Society 2006, 128 (50), 16302-16311. 27. Sakharov, D. V.; Lim, C., Zn Protein Simulations Including Charge Transfer and Local Polarization Effects. Journal of the American Chemical Society 2005, 127 (13), 4921-4929. 28. König, P. H.; Ghosh, N.; Hoffmann, M.; Elstner, M.; Tajkhorshid, E.; Frauenheim, T.; Cui, Q., Toward Theoretical Analyis of Long-Range Proton Transfer Kinetics in Biomolecular Pumps†. The Journal of Physical Chemistry A 2005, 110 (2), 548-563. 29. Khandelwal, A.; Lukacova, V.; Comez, D.; Kroll, D. M.; Raha, S.; Balaz, S., A Combination of Docking, QM/MM Methods, and MD Simulation for Binding Affinity Estimation of Metalloprotein Ligands. Journal of Medicinal Chemistry 2005, 48 (17), 5437-5447. 30. Schutz, C. N.; Warshel, A., Analyzing Free Energy Relationships for Proton Translocations in Enzymes:  Carbonic Anhydrase Revisited. The Journal of Physical Chemistry B 2004, 108 (6), 2066-2075. 31. Loferer, M. J.; Tautermann, C. S.; Loeffler, H. H.; Liedl, K. R., Influence of Backbone Conformations of Human Carbonic Anhydrase II on Carbon Dioxide Hydration:  Hydration Pathways and Binding of Bicarbonate. Journal of the American Chemical Society 2003, 125 (29), 8921-8927. 32. Cui, Q.; Karplus, M., Is a “Proton Wire” Concerted or Stepwise? A Model Study of Proton Transfer in Carbonic Anhydrase. The Journal of Physical Chemistry B 2003, 107 (4), 1071-1078. 33. Rossi, K. A.; Merz, K. M., Jr.; Smith, G. M.; Baldwin, J. J., Application of the Free Energy Perturbation Method to Human Carbonic Anhydrase II Inhibitors. Journal of Medicinal Chemistry 1995, 38 (12), 2061-2069. 34. Liang, J.-Y.; Lipscomb, W. N., Binding of substrate CO2 to the active site of human carbonic anhydrase II: a molecular dynamics study. Proceedings of the National Academy of Sciences 1990, 87 (10), 3675-3679. 35. Merz, K. M.; Banci, L., Binding of bicarbonate to human carbonic anhydrase II: A continuum of binding states. Journal of the American Chemical Society 1997, 119 (5), 863-871. 36. Piazzetta, P.; Marino, T.; Russo, N., Theoretical investigation on the restoring step of the carbonic anhydrase catalytic cycle for natural and promiscuous substrates. Arch Biochem Biophys 2015, 582, 101-106. 37. Mahinthichaichan, P.; Gennis, R. B.; Tajkhorshid, E., All the O2 consumed by Thermus thermophilus cytochrome ba3 is delivered to the active site through a long, open hydrophobic tunnel with entrances within the lipid bilayer. Biochemistry 2016. 38. Cortesciriano, I.; Bouvier, G.; Nilges, M.; Maragliano, L.; Malliavin, T. E., Temperature Accelerated Molecular Dynamics with Soft-Ratcheting Criterion Orients Enhanced Sampling by Low-Resolution Information. Journal of Chemical Theory & Computation 2015, 11 (7), 3446. 39. Regmi, C. K.; Bhandari, Y. R.; Gerstman, B. S.; Chapagain, P. P., Exploring the diffusion of molecular oxygen in the red fluorescent protein mCherry using explicit oxygen molecular dynamics simulations. The journal of physical chemistry B 2013, 117 (8), 2247-2253. 40. Shadrina, M. S.; English, A. M.; Peslherbe, G. H., Benchmarking Rapid TLES Simulations of Gas Diffusion in Proteins: Mapping O2 Migration and Escape in Myoglobin as a Case Study. Journal of 29

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

chemical theory and computation 2016. 41. Prinz, J.-H.; Wu, H.; Sarich, M.; Keller, B.; Senne, M.; Held, M.; Chodera, J. D.; Schütte, C.; Noé, F., Markov models of molecular kinetics: Generation and validation. The Journal of chemical physics 2011, 134 (17), 174105. 42. Chen, G.; Xu, W.; Lu, D.; Wu, J.; Liu, Z., Markov-state model for CO2 binding with carbonic anhydrase under confinement. The Journal of Chemical Physics 2018, 148 (3), 035101. 43. Chen, G.; Kong, X.; Lu, D.; Wu, J.; Liu, Z., Kinetics of CO 2 diffusion in human carbonic anhydrase: a study using molecular dynamics simulations and the Markov-state model. Phys. Chem. Chem. Phys. 2017, 19 (18), 11690-11697. 44. Byvatov, E.; Fechner, U.; Sadowski, J.; Schneider, G., Comparison of support vector machine and artificial neural network systems for drug/nondrug classification. Journal of chemical information and computer sciences 2003, 43 (6), 1882-1889. 45. Gini, G.; Lorenzini, M.; Benfenati, E.; Grasso, P.; Bruschi, M., Predictive carcinogenicity: a model for aromatic compounds, with nitrogen-containing substituents, based on molecular descriptors using an artificial neural network. Journal of chemical information and computer sciences 1999, 39 (6), 1076-1080. 46. Ma, X.; Li, Z.; Achenie, L. E.; Xin, H., Machine-Learning-Augmented Chemisorption Model for CO2 Electroreduction Catalyst Screening. Journal of Physical Chemistry Letters 2015, 6 (18), 3528. 47. Janet, J. P.; Chan, L.; Kulik, H. J., Accelerating Chemical Discovery with Machine Learning: Simulated Evolution of Spin Crossover Complexes with an Artificial Neural Network. The journal of physical chemistry letters 2018, 9 (5), 1064-1071. 48. Esfandiari, K.; Ghoreyshi, A. A.; Jahanshahi, M., Using Artificial Neural Network and Ideal Adsorbed Solution Theory for Predicting the CO2/CH4 Selectivities of Metal–Organic Frameworks: A Comparative Study. Industrial & Engineering Chemistry Research 2017, 56 (49), 14610-14622. 49. Puri, M.; Pathak, Y.; Sutariya, V. K.; Tipparaju, S.; Moreno, W., Artificial Neural Network for Drug Design, Delivery and Disposition. Academic Press: 2015. 50. Håkansson, K.; Carlsson, M.; Svensson, L. A.; Liljas, A., Structure of native and apo carbonic anhydrase II and structure of some of its anion-ligand complexes. Journal of molecular biology 1992, 227 (4), 1192-1204. 51. Gordon, J. C.; Myers, J. B.; Folta, T.; Shoja, V.; Heath, L. S.; Onufriev, A., H++: a server for estimating p K as and adding missing hydrogens to macromolecules. Nucleic acids research 2005, 33 (suppl_2), W368-W371. 52. Vanommeslaeghe, K.; MacKerell Jr, A. D., Automation of the CHARMM General Force Field (CGenFF) I: bond perception and atom typing. Journal of chemical information and modeling 2012, 52 (12), 3144-3154. 53. Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L., Comparison of simple potential functions for simulating liquid water. The Journal of chemical physics 1983, 79 (2), 926-935. 54. Huang, J.; MacKerell, A. D., CHARMM36 all‐atom additive protein force field: Validation based on comparison to NMR data. J Comput Chem 2013, 34 (25), 2135-2145. 55. Pronk, S.; Páll, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; Van Der Spoel, D., GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit. Bioinformatics 2013, 29 (7), 845-854. 30

ACS Paragon Plus Environment

Page 30 of 32

Page 31 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

56. Humphrey, W.; Dalke, A.; Schulten, K., VMD: visual molecular dynamics. Journal of molecular graphics 1996, 14 (1), 33-38. 57. Sjöblom, B.; Polentarutti, M.; Djinovic-Carugo, K.; Forster, R. E., Structural Study of X-ray Induced Activation of Carbonic Anhydrase. Proceedings of the National Academy of Sciences of the United States of America 2009, 106 (26), 10609. 58. Hessa, T.; Kim, H.; Bihlmaier, K.; Lundin, C.; Boekel, J.; Andersson, H.; Nilsson, I.; White, S. H.; von Heijne, G., Recognition of transmembrane helices by the endoplasmic reticulum translocon. Nature 2005, 433 (7024), 377. 59. Wang, P.-h.; Best, R. B.; Blumberger, J., Multiscale simulation reveals multiple pathways for H2 and O2 transport in a [NiFe]-hydrogenase. Journal of the American Chemical Society 2011, 133 (10), 3548-3556. 60. Berezhkovskii, A.; Hummer, G.; Szabo, A., Reactive flux and folding pathways in network models of coarse-grained protein dynamics. The Journal of chemical physics 2009, 130 (20), 05B614. 61. Wang, P.-h.; Best, R. B.; Blumberger, J., A microscopic model for gas diffusion dynamics in a [NiFe]-hydrogenase. Phys. Chem. Chem. Phys. 2011, 13 (17), 7708-7719. 62. Pedregosa, F.; Varoquaux, G.; Gramfort, A.; Michel, V.; Thirion, B.; Grisel, O.; Blondel, M.; Prettenhofer, P.; Weiss, R.; Dubourg, V., Scikit-learn: Machine learning in Python. Journal of machine learning research 2011, 12 (Oct), 2825-2830. 63. Kingma, D. P.; Ba, J., Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980 2014.

31

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

TOC graphic

32

ACS Paragon Plus Environment

Page 32 of 32