J. Phys. Chem. B 2008, 112, 1549-1554
1549
Molecular Simulations of Solute Transport in Xylose Isomerase Crystals Kourosh Malek*,†,‡ and Marc-Olivier Coppens‡,§ Physical Chemistry and Molecular Thermodynamics, DelftChemTech, Delft UniVersity of Technology, Julianalaan 136, 2628 BL Delft, The Netherlands, and Department of Chemistry, UniVersity of Tehran, P.O. Box 14155-6455, Tehran, Iran ReceiVed: December 30, 2006; In Final Form: May 30, 2007
Cross-linked enzyme crystals (CLECs) enclose an extensive regular matrix of chiral solvent-filled nanopores, via which ions and solutes travel in and out. Several cross-linked enzyme crystals have recently been used for chiral separation and as biocatalysts. We studied the dynamics of solute transport in orthorhombic D-xylose isomerase (XI) crystals by means of Brownian dynamics (BD) and molecular dynamics (MD) simulations, which show how the protein residues influence the dynamics of solute molecules in confined regions inside the lattice. In the BD simulations, coarse-grained beads represent solutes of different sizes. The diffusion of S-phenylglycine molecules inside XI crystals is investigated by long-time MD simulations. The computed diffusion coefficients within a crystal are found to be orders of magnitude lower than in bulk water. The simulation results are compared to the recent experimental studies of diffusion and reaction inside XI crystals. The insights obtained from simulations allow us to understand the nature of solute-protein interactions and transport phenomena in CLECs, which is useful for the design of novel nanoporous biocatalysts and bioseparations based on CLECs.
Introduction Cross-linked enzyme crystals (CLECs) are crystalline nanoporous materials. 1-3 The enzyme molecules are joined in a three-dimensional lattice with periodically ordered nanochannels between them, ranging in size between 0.2 and 10 nm.1 The lattice provides structural strength and durability, since the energy of crystallization as well as that of the cross-linkers is added to the energy needed to denature the enzyme.2 The nanochannels allow substrate molecules to pass through the crystal lattice and access the active sites of the enzyme molecules, where they can adsorb and possibly react. This offers great potential to use CLECs for purification in the fine chemical and pharmaceutical industry.3 The size of CLEC particles can be controlled through the details of the crystallization process. In recent years, several cross-linked protein crystals have been fabricated and used in chiral separations4,5 and biocatalysis.6,7 Understanding the nature of transport of solutes and ions in CLECs is relevant to the aforementioned biotechnological processes.1-7 In most of these applications, the size of the substrate,9 the pore network structure, and the pore size distribution11,12 of the enzyme crystals are the most important factors controlling the mass transfer inside the crystal. Because of masstransport limitations, the activity of CLECs is often reduced compared to that of soluble protein. The effective diffusion coefficients of solutes can also be considerably lower than those in bulk water.11,13,14 The diffusion of solutes within some protein crystals has been investigated using different experimental techniques13-15 and computer simulations.11,12,16 Velev et al. studied the diffusion of surfactants in lysozyme crystals using fluorescent probes by means of quantitative fluorescence microscopy.13 The measured diffusion coefficients ranged from 2 × 10-14 to †
Delft University of Technology. Tehran University. § Present address: Howard P. Isermann Department of Chemical and Biological Engineering, Rensselaer Polytechnic Institute, Troy, NY 12180. ‡
30 × 10-14 m2 s-1. By using 3D laser scanning confocal microscopy (CLSM), Cvetkovic et al. visualized the diffusion of fluorescein and rhodamine B into tetragonal, orthorhombic, and triclinic lysozyme crystals.14,15 Their results showed that the transport of solute molecules depends on the chemical nature of the solutes (e.g., hydrophobic vs hydrophilic), the charge and size of the solutes and the structure of the pore network. The diffusion of fluorescein with an average diameter of 0.69 ( 0.02 nm in orthorhombic lysozyme crystals was found to be highly anisotropic, and the diffusion coefficient along the z-axis was calculated to be around (∼7.0 ( 0.5) × 10-13 m2 s-1.15 In addition to advanced experimental techniques,13-15 molecular simulations help to correlate the reactivity of the protein and solute transport with the nanoporosity of the enzyme crystal, from the meso scale down to the atomistic level.11,12,16 In earlier work, we carried out dynamic Monte Carlo and Brownian dynamics (BD) simulations to study electrostatic and steric confinement effects on the mobility of small model spherical probes in orthorhombic and tetragonal lysozyme crystals.11 These simulations revealed a transition between the dominance of electrostatic effects on transport of small charged probes and steric confinement for larger probes. In addition to molecular simulations, differential equations describing diffusion of small molecules could be solved to relate the diffusion time to the rate of catalytic turnover in several protein crystals.14-16 Recently, Geremia et al.16 developed a simple theoretical model to reproduce the diffusion of small molecules in lysozyme crystals. The diffusion times ranged from a few seconds to several hours and were in agreement with available experimental data. A combination of multiscale molecular simulations and continuum calculations is generally needed to quantitatively study the effects of solvent and diffusing molecules. Here, we report detailed molecular simulations of diffusion of solute species in xylose isomerase crystals. D-xylose isomerase (XI) is one of the most widely used industrial enzymes, due to
10.1021/jp069047i CCC: $40.75 © 2008 American Chemical Society Published on Web 01/17/2008
1550 J. Phys. Chem. B, Vol. 112, No. 5, 2008
Malek and Coppens
its ability to catalyze the isomerization of D-glucose to Dfructose.8 Recently, cross-linked xylose isomerase (CLXI) crystals have been used to simultaneously carry out catalysis and product separation,9 by enhancing the production of fructose from glucose.10 CLXI crystals are considerably more stable than soluble XI.9,10 CLXI crystals can separate product and reactant molecules on the basis of their molecular weight, polarity, charge, and chirality.3,9,10 We apply BD simulations to study diffusion of spherical probes of different sizes inside nanopores of XI crystals. These spherical probes mimic the size of solute molecules inside water-filled channels. The anisotropy of the diffusion is also investigated. To compare diffusion of a real solute molecule to that predicted by BD simulations, independent molecular dynamics (MD) simulations were performed using S-2-phenylglycine (S-PG) as the solute. The simulation results are validated by experimental studies of diffusion and reaction inside CLXI crystals.8,9
Figure 1. Computer-generated images of an orthorhombic I222 D-xylose isomerase lattice along the (a) x, (b) y, and (c) z crystallographic axes. A single unit cell is labeled by a square. The surface representation was computed from the electron density calculations using the crystal structure data (1XIB) available in the Brookhaven Protein Data Bank (PDB). The solvent channels are shown in white. Hydrophilic, hydrophobic, and polar regions are shown in blue, red, and green, respectively.
Methodology Protein Crystal Model. The coordinate file for orthorhombic xylose isomerase crystal was obtained from the Protein Data Bank17 (http://www.rcs-b.org), with PDB file 1XIB. One XI molecule consists of 388 amino acids with 3904 atoms in total. The simulations were done at pH 7. The residues Glu and Asp were taken to be deprotonated, while Lys, Arg, and His residues were protonated.8 This leads to -20 electron charges per XI molecule. Sodium counter-ions were then added for electroneutrality. A unit cell (a ) 9.388 nm, b ) 9.968 nm, c ) 10.290 nm) of the orthorhombic crystal was first constructed from eight proteins according to the lattice symmetry I222. By translating these coordinates along the x, y, and z-axes (parallel to the unit cell crystallographic directions), a box of unit cells was formed, as shown in Figure 1a. In our MD simulation, the coordinates of 378 crystallographic water molecules per protein were added to the unit cell as well. BD Simulations. The effective diffusion coefficient of spherical, neutral probes within the crystal channels is calculated by BD simulations.11,18 The initial simulation box had a size of 2 × 2 × 2 unit cells (a ) 18.776 nm, b ) 19.936 nm, c ) -20.580 nm), and was increased by using periodic boundary conditions. In our BD simulations, we removed the explicit solvent molecules and replaced them with an implicit continuum solvent medium. We also ignored internal motions of the XI molecules, and therefore, the XI crystals were modeled as rigid bodies. This ignores the effect of protein fluctuations on transport inside the protein channels, which could matter, especially for solutes of sizes close to the pore diameter. Previous MD simulations showed that the protein atoms only fluctuate slightly, and that these fluctuations are unlikely to cause significant changes in pore size.12 However, this may only be true for the transport of small probes (e.g., glucose) of sizes considerably smaller than the channel diameter. Another drawback to BD simulations is the lack of single-point contacts
(e.g., hydrogen bonding and dipole-dipole interactions) between the real solute molecules and the protein side chains. Such interactions are ignored, where using coarse-grained spherical solutes. Nonetheless, our BD simulations provide insight into size-exclusion and steric effects on the diffusion of solutes in XI crystals. The latter is computationally expensive and difficult to achieve for large solute molecules by means of full-atom molecular dynamics simulations. A single polar, spherical bead represents the solute in our BD simulations. The non-bonded interactions between solute beads and protein atoms were modeled by a Lennard-Jones (LJ) potential:
[( ) ( ) ]
ULJ(r) ) 4ij
σij r
12
-
σij r
6
(1)
The effective solute diameter (σij) was varied from 0.2 to 2.2 nm. The LJ interactions C6 (4σ6) and C12 (4σ12) for beads with different sizes σ were taken from the strength of nonbonding interactions in the GROMACS database, which are illustrated in Table 1. One solute bead was simulated at a time during each run, and, therefore, solute-solute interactions were not taken into account in our BD simulations. The probe-probe interaction is given in Table 1 for the sake of comparison. The short-range electrostatic interactions (counterion-counterion and counterion-residues) in our simulations were calculated using the Coulombic potential. Multiple additional runs with different starting configurations were performed. The motion of the ith probe with mass mi and charge qi was governed by the Langevin equation
TABLE 1: Nonbonding Interaction Parameters in Terms of C6 (4Eσ6) and C12 (4Eσ12) for Probes with Different Sizesa P-P
P-C
P-O
P-N
probe size (nm)
C6
C12
C6
C12
C6
C12
C6
C12
0.2 0.6 1 1.4 1.8 2.2
0.00128 0.93312 20 150.5907 680.2445 2267.598
8.19E-08 0.043536 20 1133.878 23136.63 257100.1
0.000289 0.023448 0.288581 1.684808 6.575819 19.93621
1.16E-08 7.64E-05 0.011567 0.394247 6.00575 55.20173
0.000404 0.03843 0.496677 2.969677 11.75541 35.97351
1.20141E-08 0.000108591 0.018138807 0.648454582 10.16100629 95.1539386
0.000525 0.04862 0.623191 3.711177 14.65581 44.77904
1.64E-08 0.000141 0.023117 0.819812 12.78529 119.3549
a
Key: P, probe; C, carbon; O, oxygen; N, nitrogen.
Solute Transport in Xylose Isomerase Crystals
mi
dVi ) -miγiνi + FR(t) + qiζi dt
J. Phys. Chem. B, Vol. 112, No. 5, 2008 1551
(2)
The first term on the right-hand side of eq 2 corresponds to an average frictional force with a friction coefficient miγi. FR(t) represents the random part of the collisions and rapidly fluctuates around a zero mean, while qiζi is the electrical force acting on an ion. The properties of the system were extracted from time averages taken over atomistic trajectories, which are stochastic in the case of BD.18 In our Brownian dynamics simulations, long-range electrostatic interactions were calculated using the particle mesh Ewald (PME) method with a grid spacing of 0.12 nm and fourth-order interpolation.19 The medium outside a cutoff sphere (radius RC) centered at a charge group was modeled as a dielectric continuum of specified relative permittivity equal to that of the bulk solvent, r ) 78. Brownian motion of the solute probes was simulated using the GROMACS package20,21 with a time step of 0.25 ps. All simulations were performed under constantvolume. MD Simulations. We used MD simulations to examine diffusion of S-PG in a fully hydrated orthorhombic XI crystal. Simulations were carried out using a fully atomistic 1 × 1 × 2 XI lattice (a ) 9.388 nm, b ) 9.968 nm, c ) 20.580 nm), Figure 2. The lattice model contained S-PG molecules (Scheme 1), sodium counterions, and water. The total system consisted of 16 XI molecules, 7 S-PGs, 6048 crystallographic water molecules, 32946 noncrystallographic water molecules and 320 sodium ions, leading to a grand total of 279475 atoms. First, the system was equilibrated for 500 ps using harmonic position restraints (1000 kJ mol-1 nm-2). Then the production runs were performed for another 7.5 ns, the last 7 ns of which were used in the analysis. Simulations were carried out using the GROMOS96 force field.22,23 In our force field, interactions between atoms are divided into nonbonded interactions, between any pair of atoms that are within a given cutoff radius, and bonded interactions between atoms connected by chemical bonds. In the case of the nonbonded interactions (electrostatic and van der Waals), a partial charge and parameters for repulsion and attraction are assigned to each atom. Bonded interactions consist of bond, angle, and dihedral terms. Bond bending and angle bending are given by simple harmonic potentials. The torsionalrotational potential for the dihedral angle is a periodic function with a 3-fold barrier. The partial charges and interaction parameters for all species were taken from the GROMOS96 database.23 In our MD simulation, a cutoff of 1.4 nm was used for van der Waals interactions, while a cutoff of 1.0 nm and PME with a grid spacing of 0.12 nm and fourth-order interpolation were used for electrostatic interactions.19 During the simulation, the potential energy and the total energy were monitored to check whether the system was in equilibrium. The root-mean-square deviation of the atoms from the initial configuration was used to determine the equilibrium and stability of the system. MD simulations were performed in a canonical (NVT) ensemble at 300 K. Similar to BD simulations, the temperature was controlled by the Berendsen weak-coupling algorithm, separately for protein, S-PG, and water plus ions with a time constant of 0.1 ps and a temperature of 300 K. During the simulations, polar hydrogen atoms were treated as dummy atoms with an increased mass of 4 Da. This allowed the integration time step to be 5 fs, which improves the simulation efficiency.24,25 It should be noticed, however, that the redistribution of atomic masses increases the momentum of polar hydrogen atoms in the protein.
Figure 2. View of the simulation box along the z axis (a) and y axis (b).
SCHEME 1
As a result, the rotational motion is slowed down, which may have considerable impact on translational diffusion as well. The diffusion coefficient appears to be reduced by a factor of 10.24 A smaller time step of 1-2 fs is a reasonable choice, which engages a higher degree of freedom for S-PG molecules. We are interested in effective dynamic properties of S-PG molecules based on center-of-mass diffusion. The choice of a larger time step allows one to achieve longer simulation times, and, therefore, we consider the modified polar hydrogen and the 5 fs time step to be an acceptable approximation. Simulations were carried out with the GROMACS package20,21 (http://www.gromacs.org). Results were visualized using the VMD v1.8.326 commercial package. Results and Discussion Figure 1 shows the atomic model of a 2 × 2 × 2 unit cell of the XI crystal structure. The main, widest pore (PZ) lies along the z axis, Figure 1a. There are other pores perpendicular to the PZ pore, namely PX and PY, which run along the x- and y-axes, respectively, Figures 1b,c. The active sites of the D-XI enzyme for the isomerization of D-glucose are usually defined by residues His54, Phe94, Lys183, Asp287, Glu181, and His220.27 It should be noticed that not all of these residues are accessible from within the pore space in a CLXI crystal. Some of the residues lie along the wall of the pores, while many of them are buried inside the lattice. Accordingly, the activity of the XI crystals for the isomerization reaction of D-fructose depends on the accessibility of these residues inside the channels. However, all pores can be considered as diffusion channels. The pore network within the XI lattice is determined by using
1552 J. Phys. Chem. B, Vol. 112, No. 5, 2008
Malek and Coppens
Figure 3. Pore structure and pore radius profile along the pore axis for PZ and PX pores (black, PZ; red, PX). Hydrophilic and hydrophobic regions are shown in blue and in red, respectively.
Figure 4. (a) Simulated diffusivities of an uncharged bead in an orthorhombic I222 D-xylose isomerase crystal as a function of the sphere radius. (b) Diffusion anisotropy: ratio of the diffusivity along the z axis to the three-dimensional diffusivity as a function of the probe radius.
a procedure explained elsewhere,11 based on the HOLE28 and CHANNEL29 algorithms. The radius of a pore is determined at any given distance along the pore axis by calculating the maximum size for a spherical probe.29 Figure 3 illustrates the pore shape, along with the pore radius profile for the PZ and PX pores. The pore profile along the y-axis (PY) is very similar to that along the x-axis (PX). PZ pores form straight channels and contain active sites Asp287 and His54 on the wall. The profile of the pore size along the pore axis shows that there are constricted regions inside the channel. The pore radius in PZ slowly decreases from over 2.5 nm to slightly less than 1.4 nm at its narrowest part. For PX and PY, on the other hand, pores can be represented as consisting of big cavities with narrow cylindrical interconnections. The pore radius in PX varies from 0.9 to 2.0 nm and from 1.0 to 1.6 nm for PY. Figure 4a shows diffusion results from BD simulations for noncharged solutes of radius 0.2-2.2 nm in the orthorhombic D-XI crystal. Clearly, the diffusivity drops sharply when the radius of the mobile bodies approaches the channel size. Note that probes with a radius of up to about 2 nm are still able to diffuse through the pore space, even if the earlier Monte-Carlo generated pore radius profiles in Figure 3 seem to suggest that this should be impossible. The reason for this is that the pores do not have a circular cross-section, and they do not intersect in spherical intersections. They form a three-dimensional
Figure 5. (a) Displacement of seven S-PG molecules in the pore regions of the XI crystal along different unit cell coordinates. Different colors represent different S-PG molecules. (b) MSD of S-PG molecules in the XI crystal versus time. The inset indicates the size of the diffusing molecules compared to the channel diameter.
network, and have a convoluted surface with cavities along the walls. Brownian Dynamics simulations therefore reflect more accurately the real accessibility of this pore space to a diffusing probe, and Figure 4 is a good measure for this accessibility. Molecules larger than 2 nm in radius cannot easily enter pores because their size is very close to the crystal pore size. The flexibility of the protein might play an increasingly important role in the diffusion of probes larger than about 1.5 nm in radius. For larger radii, the 3D matrix of crystal channels is reduced to a quasi-1D system of isolated tubes along the z direction. These calculations only account for steric limitations of the crystal to diffusion, but other factors might influence diffusion in real cross-linked crystals as well. The presence of cross-linker agent within the channels of the real crystal is one of the additional obstacles.1,12 The most important source reducing the mobility of solutes is their binding to surface groups of protein molecules. Such interactions can be better described in the course of detailed, fully atomistic calculations, similar to our MD simulation for diffusion of S-PG in XI crystals. Diffusion in the crystals may also be significantly anisotropic. The ratio of the diffusivity along the z axis to the three-dimensional diffusivity, D/D3D, strongly depends on the size of the probe, Figure 4b. For probes of sizes > 1.25 nm, a transition from three-dimensional diffusive motion to one-dimensional could be responsible for higher D/D3D values than one. The total diffusivity and the diffusivity in the z direction are not identical for smaller probes, because there are channels that contribute significantly to diffusion in other directions, Figure 1. The diffusion anisotropy is most pronounced for larger probes whose size is between 1 and 1.4
Solute Transport in Xylose Isomerase Crystals
J. Phys. Chem. B, Vol. 112, No. 5, 2008 1553
nm, the diffusivity along the z axis and the three-dimensional diffusivity coincide for these large probes, since only transport along the z-axis is possible. To study diffusion of a solute (S-PG) inside the XI crystal pore space, we have performed MD simulations of S-PG transport in a fully hydrated XI crystal. Protein atoms are allowed to fluctuate. The solute concentration was kept at 3M, which is in the range of the typical sugar concentration used in previous experiments.5-10 This concentration corresponds to seven solute molecules per simulation box. Both our simulations and previous experimental studies have shown that diffusion in protein crystals depends on solute concentration.16 Our assumption that seven diffusing molecules can interpret the MSD behavior of the S-PG molecules is obviously a gross simplification. We believe, however, that the present system size is sufficient to understand the principle mechanism of solute diffusion and solute-protein interactions in XI crystal. We also randomly incorporated the position of the seven S-PG molecules inside the PZ pore, Figure 2a. The dynamical motion of the solute molecules was sampled every 1 ps during the 7 ns simulation. Figure 5a shows the displacement of the seven S-PGs along the x, y, and z-axes of the simulation box, as a function of time. The location of the solute molecules is sufficiently randomized during the simulation time, so that the diffusion properties should not significantly depend on the arbitrary initial placement of the S-PG molecules. The motion of S-PG molecules occurs in many jumps, with little motion between the jumps. A few S-PG molecules travel within the pore space, with displacements on the order of 10 nm, while a few other ones remain around the same axial (z) position in the pore region, moving only slightly in the xy-plane. Overall, S-PG molecules move freely; some travel all the way within the pore network, a few remain in the pore, and some go deep into the pore and return after some time. It was also observed that at some points along the trajectory, a group of two or three S-PG molecules establish hydrogen bonding with the protein residues or perform a collective motion as a united body. The latter was especially observed at the longer simulation times, which creates an artifact in the diffusion analysis of single S-PG molecules. To calculate the self-diffusion D of S-PG molecules in the pore, the Einstein relationship is used:
〈|∆rj| 〉t0 2
D ) lim tf∞
6t
(3)
where 〈|∆rj|2〉t0 is the mean square displacement (MSD) of S-PG molecules during a time t, averaged over the ensemble of molecules in the pore, from an initial moment, t0 on that their dynamical motion is followed in the pore.30 In MD simulations, t has to be large compared to the relaxation times of microscopic processes, which corresponds to the correlation time of the velocity autocorrelation function.31,32 For the determination of the S-PG diffusivity in our MD simulations in the presence of water and ions, each time increment ∆t was assumed to be 50 ps. At the start of each new time frame of 50 ps, the contribution of each S-PG molecule to the mean square displacement over the current time frame was calculated. Afterward, the squared displacement was averaged over all molecules and time origins. The MSD at each time interval was then plotted over simulation time (7 ns). If the displacement reflects a purely Brownian motion, the choice of time frames is arbitrary. In reality, however, it must be sufficiently large to incorporate all dynamical memory effects present in solution.33 Moreover, a linear dependence of the MSD with t at large times has to be verified. Figure 5b shows the MSD vs time for S-PG molecules
diffusing in the pore network of the XI crystal. Only trajectories of individual S-PG molecules are included in calculation of the MSD. The collective motion of S-PG molecules is not included in Figure 5b. Overall, the log-log behavior of the MSD vs time is linear with a slope close to 1. This shows that in a fully hydrated pore of the XI lattice, the Einstein relation, eq 3, can describe the diffusion of S-PG molecules. Let us approximate each S-PG molecule by a spherical probe with a radius of 0.25 nm. The calculated diffusion coefficient of S-PG in pore PZ (2.3 × 10-11 m2 s-1) from MD simulations is comparable to the diffusivity of a probe of size 0.25 nm, obtained from BD simulations (1.88 × 10-11 m2 s-1). The deviation is most likely a result of single-point interactions between the solute and side chains of the XI, which are usually characterized by hydrogen bonds and dipole-dipole interactions. Our preliminary analysis based on hydrogen bonds, the number of contacts, and the minimum distance between S-PG molecules and residues on the pore wall showed that the S-PG molecules spend considerable time in the vicinity of XI active site residues His54, Asp287, and Lys183. The actual activity of CLXI crystals for the isomerization of D-glucose to D-fructose has recently been measured by Vuolanto and Leisola.34 It was shown that the activity decreased by increasing the crystal size, and the effective diffusivity of the isomerization reaction, De, was estimated from the experimental results.32 The Thiele modulus φ is a dimensionless number, representative for the ratio of the intrinsic chemical reaction rate in the absence of mass transfer limitations to the rate of diffusion through a catalyst particle. The Thiele modulus is calculated from the internal effectiveness factor η for spherical catalyst pellets, by assuming a first-order isomerization reaction:
η)
3 [Φ coth(Φ) - 1] Φ
(4)
The Thiele modulus is given by
Φ)R
x
rF DeC
(5)
in which De is the diffusion coefficient, R is the crystal size, F is the crystal density, C is the reactant concentration, and r is the reaction rate. The diffusion coefficient thus obtained ranged from 0.5 × 10-11 to 1.2 × 10-11 m2 s-1 by fitting the reactivity data to the diffusion equation.34 Our BD simulations for a probe of a size comparable to that of the hydrated D-glucose (0.42 nm) yield a diffusivity of around 1.55 × 10-11 m2 s-1, close to the experimental values.34,35 The calculated effective diffusivity of D-glucose from our simulations is more than one order of magnitude lower than the diffusivity of D-glucose or D-fructose in water at 25 °C (∼5.2 × 10-10 m2 s-1).36 For glycine, the experimental value of the diffusivity in water at 300 K (1.21 × 10-9 m2 s-1) is more than 2 orders of magnitude higher than the diffusivity in XI crystals.35,36 Note that, in protein crystals, diffusion may be affected by the charge on the protein (which changes with pH), the polarity of the probe, and the size of the molecules and pore channels. In general, the intracrystalline diffusivities were found to be about 1-2 orders of magnitude lower than those of the corresponding free molecules in water.34-36 The diffusion of solutes (radius up to 0.6 nm) within orthorhombic lysozyme crystals, calculated from our previous BD simulations [(0.5-1.5) × 10-10 m2 s-1], appears to be slightly faster than that in XI crystals. Despite smaller pore sizes in lysozyme crystals, and in addition to the wetting properties of the lattice, the deviation is most likely a result of the highly
1554 J. Phys. Chem. B, Vol. 112, No. 5, 2008 oriented pore network structure of the orthorhombic lysozyme crystals. The latter leads to more anisotropic diffusion inside lysozyme, compared to the motion in the three-dimensional pore network in XI crystals, which is more isotropic for solutes in that size range. Both in our simulations and in experiments, diffusion coefficients are charge dependent and decrease with increasing molecular weight (volume). Although our simulations only dealt with the effect of the crystal packing and charge density on diffusion, they could be used to model chiral separation processes in XI crystals. Chromatographic separation is often based on the interaction of molecules in a mixture with the active centers of the enzyme molecules in the lattice and could therefore be simulated in a manner similar to the one proposed here. Conclusions In the present study, we performed independent BD and MD simulations of solute diffusion in nanopores of orthorhombic XI crystals. XI crystals consist of a three-dimensional solventfilled pore network with pore radii between 0.9 and 2.5 nm. BD simulations of diffusion of coarse-grained beads of radii between 0.2 and 2.2 nm through these nanopores provided insight into electrostatic and steric confinement effects on the solute transport in XI crystals. To complement the results from BD simulations, the stereochemical interactions between protein residues and a real solute, S-PG, were investigated by means of fully atomistic MD simulations. Our analysis showed that the S-PG molecules mostly interact with residues His54, Asp287, and Lys183. In general, the diffusivities of solute species were found to be 1 to 2 orders of magnitude lower than those of the corresponding free molecules in water. These simulations suggest that steric (size exclusion) and electrostatic (dipoledipole and hydrogen bonding) interactions between dipolar solutes and amino acids exposed on the pore wall are mostly responsible for this reduction. Acknowledgment. Financial support from the LifeTech program of Delft University of Technology is gratefully acknowledged. We thank Dr. Antti Vuolanto and Prof. M. Leisola for useful discussions. References and Notes (1) Margolin, A. L.; Navia, M. A. Angew. Chem., Int. Ed. 2001, 40, 2205. (2) Vilenchik, L. Z.; Griffith, J. P.; Clair, N. St.; Navia, M. A.; Margolin, A. L. J. Am. Chem. Soc. 1998, 120, 4290. (3) Margolin, A. L.; Vilenchik, L. Z. U.S. Patent WO9813119, 1996. (4) Margolin, A. L.; Persichetti, R. A.; Clair, N. L. St.; Khalaf, N. K. U.S. Patent 6140475, 2000. (5) Vuolanto, A.; Leisola, M.; Jokela, J. Biotechnol Prog. 2004, 20, 771.
Malek and Coppens (6) Vuolanto, A.; Pastinen, O.; Schoemaker, H. E.; Leisola, M. Biocatal. Biotransform. 2002, 20, 235. (7) Vuolanto, A.; Pastinen, O.; Schoemaker, H. E.; Leisola, M. Biocatal. Biotransform. 2002, 20, 235. (8) Vangrysperre, W.; Christophe, A.; Kersters-Hillderson, H.; Tempst, P. Biochem. J. 1989, 263, 195. (9) Leisola, M.; Jokela, J.; Finell, J.; Pastinen, O. Biotechnol. Bioeng. 2001, 72, 501. (10) Jokela, J.; Pastinen, O.; Leisola, M. Enzyme Microb. Technol. 2002, 31, 67. (11) Malek, K.; Odijk, T.; Coppens, M.-O. ChemPhysChem 2004, 5, 1596. (12) Malek, K.; Odijk, T.; Coppens, M.-O. Nanotechnology 2005, 16 S522. (13) Velev, O. D.; Kaler, E. W.; Lenhoff, A. M. J. Phys. Chem. B 2000, 104, 9267. (14) Cvetkovic, A.; Picioreanu, C.; Straathof, A. J. J.; Krishna, R.; van der Wielen, L. A. M. J. Am. Chem. Soc. 2005, 127, 875. (15) Cvetkovic, A.; Straathof, A. J. J.; Hanlon, D. N.; van der Zwaag, S.; Krishna, R.; van der Wielen, L. A. M. Biotechnol. Bioeng. 2004, 86, 389. (16) Geremia, S.; Campagnolo, M.; Demitri, N.; Johnson, L. N. Structure 2006, 14, 393. (17) Carrell, H. L.; Hoier, H.; Glusker, J. P. Acta Crystallogr., D: Biol. Crystallogr. 1994, 50, 113. (18) Ermak, D. L.; McCammon, J. A. J. Chem. Phys. 1978, 69, 1352. (19) de Leeuw, S. W.; Perram, J. W.; Smith, E. R. Proc. R. Soc. London, Ser. A 1980, 373, 27. (20) Lindahl, E.; Hess, B.; van der Spoel, D. J. Mol. Model. 2001, 7, 306. (21) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. Comput. Phys. Commun. 1995, 91, 43. (22) van Gunsteren, W. F.; Kruger, P.; Billeter, S. R.; Mark, A. E.; Eising, A. A.; Scott, W. R. P.; Heneberger, P. H.; Tironi, I. G. The GROMOS96 Manual and User Guide, Biomos and Hochschulverlag AG an der ETH Zu¨rich: Groningen, The Netherlands 1996. (23) van Gunsteren, W. F.; Berendsen, J. C. Angew. Chem., Int. Ed. Engl. 1990, 29, 992. (24) Feenstra, K. A.; Hess, B.; Berendsen, H. J. C. J. Chem. Phys. 1999, 20, 786. (25) van der Spoel, D.; van Maaren, P. J.; Berendsen, H. J. C. J. Chem. Phys. 1998, 108, 10220. (26) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graphics 1996, 14, 33. (27) Whitaker, R. D.; Cho, Y.; Cha, J.; Carrell, H. L.; Glusker, J. P.; Karplus, P. A.; Batt, C. A. J. Biol. Chem. 1995, 270, 22895. (28) Smart, O. S.; Neduvelil, J. G.; Wang, X.; Wallace, B. A.; Sansom, M. S. P. J. Mol. Graphics 1996, 14, 354. (29) Kisljuk, O. S.; Kachalova, G. S.; Lanina, N. P. J. Mol. Graphics 1994, 12, 305. (30) Bizzari, A. R.; Cannistraro, S. Phys. ReV. E 1996, 53, R3040. (31) McQuarrie, D. A. Statistical Mechanics; Harper and Row: New York, 1976. (32) Ahlstrom, P.; Teleman, O.; Jonsson, B. J. Am. Chem. Soc.1988, 110, 4198. (33) Roux, B. Biophys. J. 1998, 74, 2744. (34) Vuolanto, A. Cross-linked Protein Crystal Technology in Bioseparation and Biocatalysis Applications. Ph.D. Thesis, Helsinki University of Technology, 2004. (35) Umecky, T.; Kuga, T.; Funazukuri, T. J. Chem. Eng. Data 2006, 51, 1705. (36) van de Ven-Lucassen, I. M. J. J.; Kerkhof, P. J. A. M. J. Chem. Eng. Data 1999, 44, 93.
