International Conference on Gas Hydrates (ICGH) (6th : 2008)


Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


59278-5510.pdf [ 374.67kB ]
JSON: 59278-1.0041134.json
JSON-LD: 59278-1.0041134-ld.json
RDF/XML (Pretty): 59278-1.0041134-rdf.xml
RDF/JSON: 59278-1.0041134-rdf.json
Turtle: 59278-1.0041134-turtle.txt
N-Triples: 59278-1.0041134-rdf-ntriples.txt
Original Record: 59278-1.0041134-source.json
Full Text

Full Text

SIMULATING THREE-DIMENSIONAL GAS HYDRATE GROWTH AND INHIBITION Brent Wathen and Zongchao Jia Department of Biochemistry Queen's University Kingston, CANADA Virginia K. Walker Department of Biology Queen's University Kingston, CANADA ABSTRACT The economic and safety hazards associated with the ability of gas hydrates to form in pipelines have prompted our interest in the inhibition of hydrate growth. Antifreeze proteins (AFPs) adsorb to ice surfaces and certain AFPs can also inhibit the growth of hydrates formed from water molecules organized in cage-like formations around a central small gas molecule. A Monte Carlo computational method for simulating the growth of ice crystals has been developed and it has proved useful in the understanding of the inhibition mechanism of these proteins. We have modified this crystal growth software in order to simulate the growth of large structure II gas hydrates, consisting of millions of water and gas molecules. This represents a first step towards investigating the effectiveness of novel compounds to inhibit hydrate growth in silico. Here, we describe these software modifications, and our efforts to incorporate type I AFP molecules into the hydrate growth simulations. Because both the docking interaction and inhibition mechanism for AFP towards hydrates remains unknown, we have set up a number of inhibitor screens to investigate possible AFP-hydrate docking models. Our goal is to reproduce the changes to gas hydrate morphology that have been observed in the presence of AFP, which will guide our choices for the binding alignment between AFPs and hydrates. This alignment will be instrumental for determining the AFPI-inhibition mechanism and should prove invaluable for the development of novel, hyperactive hydrate inhibitors. Keywords: gas hydrates, antifreeze proteins, kinetic inhibitors, Monte Carlo simulations  INTRODUCTION Gas hydrate formation at well heads and in pipelines impedes the flow of oil and gas [1], increases the cost of flow assurance [1-3] and adds significant hazard to the process of oil recovery [2]. These icelike crystalline structures encapsulate gas molecules in water cages. Three distinct forms of gas hydrates are known to occur naturally [4], including the cubic structures I and II (sI and sII), and the much rarer hexagonal structure sH. All three are built from a combination of basic 512 "small" cage building blocks (12-faced pentagonal dodecahedron  structures) and other water cage structures. In the case of the sII structure, 16 of the 512 small cages combine with 8 51264 "large" cages in a periodic manner to produce a crystalline structure that has a unit cell with 136 water molecules and 24 small and large gas-trapping cages (Figure 1A). While it is not possible to predict with certainty which of these three structures will form if multiple gases are present [5], in general smaller gas molecules tend to promote sII formation, while larger ones tend to promote sH formation [6]. Although gas hydrates are non-stoichiometric compounds, the gas  Figure 1. Simulation elements. A) Hydrate sII structure, showing two gas molecules (grey spheres) trapped in large water cages, with small water cages nearby. B) Type I AFP. The alanine residues on the binding face are shown in blue. occupancy rate has been measured approaching 100% [7]. To decrease hydrate formation, two different classes of hydrate inhibitors have been employed. The first class includes widely used thermodynamic inhibitors such as alcohols and glycols to minimize hydrate growth and enhance petroleum flow [1]. Large quantities of these inhibitors are required, however, making the use of these types of inhibitors costly [1]. A second class of inhibitors, low dosage hydrate inhibitors (LDHIs), appears to disrupting gas hydrate growth following formation [8]. Antifreeze proteins (AFPs), which have long been studied because of their ability to inhibit ice growth in certain cold tolerant organisms, also act as LDHIs [9]. In retrospect, this observation is surprising given the differences in ice and hydrate structures. A detailed understanding of gas hydrate nucleation and growth would greatly assist in the development of strategies to better manage hydrate formation and perhaps uncover the mechanism(s) of LDHI and AFP hydrate growth inhibition. Molecular dynamics (MD) simulations of hydrate formation have shown that structures resembling hydrate sII can be present at early stages of hydrate sI formation [10], confirming experimental reports. Several MD studies have examined gas hydrate stability [11, 12], while others have investigated the mechanisms of LDHI inhibition [1, 8, 13]. These latter simulations have shown that some LDHIs function by binding to gas hydrates and arresting  subsequent growth [13], while others can affect hydrate growth from a distance, organizing surrounding water structure in a way that is inconsistent with gas hydrate formation [8]. Because of limited computational resources, these simulations are necessarily restricted to small numbers of molecules, perhaps hundreds or thousands of water and gas molecules, and one, or at most a few, inhibitor molecules. In particular, they are not suitable for investigating the inhibitory effects from hundreds or thousands of inhibitors working in cooperation. Previously, we have investigated the growth of hexagonal Ih ice using a novel Monte Carlo (MC) computational approach that considers the growth of whole crystals rather than single unit cells or isolated surfaces [14, 15]. This approach redistributes computational resources, effectively trading off mathematical detail for scale: the precision of MD simulations with regards to atomic positions and interactions is simplified by the use of probability, with the computational gains used to investigate ice crystals containing millions of water molecules undergoing billions of individual water association and dissociation events. This model was further expanded to allow for the inclusion of inhibitor molecules, embodied with realistic (though fixed) shape and interaction energetics. Using this model, AFP ice-growth inhibition mechanisms were investigated, and a novel inhibition theory was posited based on AFP-ice binding orientation, rather than AFP-ice binding strength [14, 15]. We have now adapted our MC technique to the study of gas hydrate formation. Here we report the progress on this project to date, including our efforts to parallelize our algorithm, and our initial efforts to identify successful AFP-hydrate binding orientations. METHODS Whole-Crystal Monte Carlo Simulations The whole-crystal MC simulation algorithm is described in refs [14, 15]. Briefly, however, the underlying model for the technique is based on the Kinetic Ising model as described by Gilmer [16]. Central to our approach is a simulation space (Figure 2A), which is a lattice of possible molecular locations and a description of connectivity between these locations. Two positions that are connected are  said to be neighbours. The structure of the simulation space is problem-specific and is chosen to reflect the molecular structure of a crystal under investigation. The simulation begins by "placing" molecules into some initial subset of the simulation space locations, with the proviso that any two of these molecules must be connected to one another via a path of connected, occupied simulation space locations. In this way, the simulation space will contain a single "crystal" (ie. the structure made up of all occupied simulation locations). Each step of the simulation will make one of two possible changes to the simulations space: either a molecule on the crystal surface is removed (a dissociation event), or a molecule is added to an interface location adjacent to a surface molecule (an association event). Surface locations are those locations in simulation space that are occupied and connected to unoccupied locations, while interface locations are those that are unoccupied and connected to occupied locations. Probability is used to determine what type of event takes place in each step. Association events happen with probability Passoc: Passoc = Locinterface * Pon  (1)  where Locinterface is the total number of interface locations, and Pon is a fixed probability of a molecule joining the crystal. Dissociation events happen with probability Pdissoc: Pdissoc = Poff ( neighbour occupancy )  (2)  where Poff, the probability of dissociation at a particular location, is a function of its neighbour occupancy. Poff ( 0 ) is set to 1.0, ensuring that any molecule that becomes "stranded" is removed from the simulation, while Poff ( full occupancy ) is set to 0.0, ensuring that "buried" molecules in crystal interiors do not dissociate. Poff values for intermediary neighbouring occupancies reflect the number of bonds that need to be broken: Poff ( i ) = exp[ -i φ / kB T ]  (3)  where i is the neighbour occupancy (ie. the number of bonds that need to be broken), φ is the bond  Figure 2. MC Simulations. A) A simplified two-dimensional simulation space, with leftright-up-down connectivity. Occupied (blue circle), surface (yellow), buried (cyan) and interface (green) positions are shown. B) A "complex" simulation molecule consists of its volume set (left, yellow squares), and its energetics set (right, depicting the interactions with each of its neighbours, determined individually). C) During simulations, the interaction energy between a "complex" molecule and crystal is the sum of the predetermined interaction energies between energy, kB and is its Boltzmanns's constant and T is molecule occupied neighbours. temperature. It is this incorporation of connectivity into the probability of dissociation that enables the specifics of a particular crystal structure to affect morphology and other growth characteristics during simulations. The basic model has been extended to allow for the incorporation of inhibitor molecules into simulation spaces during simulations. Inhibitor molecules differ from crystal molecules in two critical aspects. First, they are much larger than crystal molecules, typically occupying hundreds of times more volume than crystal molecules. Secondly, their energetics are far more complex, with non-uniform attractive and repulsive forces between themselves and neighbouring crystal molecules across their surfaces. To maintain the spatial and energetic qualities of inhibitors, a new type of simulation molecule, termed a complex molecule, was introduced into the model. Complex molecules are pre-built from a docking between a molecule of interest (ie. the inhibitor) and the crystal, from which a set of specific, connected simulation spaces representing the volume of the complex molecule is predetermined (Figure 2). Moreover, the energetics between a docked complex molecule and all possible neighbouring simulation space locations  are also predetermined (Figure 2). Inhibitors are thus incorporated into simulations as crystal molecules with two additions, (1) a volume set describing the inhibitor’s shape in the bound conformation and (2) an energetics set describing its interactions with adjacent simulation shape positions. Again, there is a fixed association probability Passoc (reflecting inhibitor concentration, diffusion through the liquid medium and the chance of a correctly oriented collision with the crystal) and a dissociation probability Pdissoc reflecting the specific energetics between a bound inhibitor molecule and those of its simulation space neighbours that are presently occupied. The Pdissoc values typically reflect the bonding energy as given in eqn. (3), while the Passoc value is generally set to a suitably low value chosen by empirical observations of simulation results. It must be emphasized that a complex molecule is static in simulations: once bound to a particular crystal plane, it cannot move nor change its structure relative to the underlying crystal. Our model efficiently maintains both the complete three-dimensional structure of a crystal in simulation space and those simulation space locations that comprise the set of crystal surface and neighbouring interface positions. In consequence, simulations containing millions of crystal molecules, and hundreds or thousands of inhibitors, can be undertaken spanning billions of individual association and dissociation events. With simulations on this scale, the cooperative effect of large numbers of inhibitor molecules on crystals, some with diameters on the order of 0.1-0.2 microns, can finally be explored. Model Implementation Changes Required for Gas Hydrate Simulations In order to adapt the model for studying gas hydrate growth, a number of changes have been made. Generally, these are the result of (i) differences in the underlying crystal structures of ice and hydrate, and (ii) the lack of a detailed LDHI-hydrate docking model. Hydrate sII has a structure that is inherently more complex than that of ice Ih due to the presence of the trapped gas molecules. Because of the fundamental role of the gas in hydrate formation, the algorithm had to be modified for multiple types of  crystal molecules, each given an independent association parameter. This model enhancement allows an investigation of the effects of gas concentration on crystal growth rates and gas occupancy rates, and thus provided a simple test of model correctness: gas hydrates crystals should melt in simulations without gas molecules. More generally, the larger, more complex sII unit cell required a redesign of the data structures used to describe the crystal, and of those used to manage the simulation space. A more significant obstacle to studying hydrate formation in the presence of inhibitors is the lack of published experimental evidence supporting any specific docking model between a gas hydrate crystal and an inhibitor. Without such a model, meaningfully oriented inhibitor molecules cannot be created for inclusion in MC simulations. To circumvent this problem, we attempted a large-scale screening process of possible docking models between the inhibitor under study, Type I AFP (AFP-I), and hydrate sII, in order to identify factors that promote inhibitor activity. Even though only a small fraction of possible docking models are being investigated, the screening process remains computationally demanding. For this process to be computationally tractable, we have had to redesign the algorithm to operate across multiple processors. The overall methodology of the MC simulation model does not lend itself well to temporal parallelization, primarily because each simulation step is dependent on the outcome of all prior simulations steps, necessitating a serial treatment. Rather than temporal, our approach to distributing simulation workload across multiple processors is spatial: we partitioned the simulation space into symmetrical regions, and assign one processor to manage the simulation in each region. Two complications from this approach arise. First, because association and dissociation events can affect the status of neighbouring locations, there are complications at regional boundaries. To address this issue, we implemented a message passing scheme between processors to ensure that all changes to the simulation space are properly coordinated across the regions. The second complication that arises from our spatial partitioning of the simulation space pertains to complex molecules: a bound inhibitor may residue in multiple  regions because of their extended volume. The association or dissociation of an inhibitor would thus be a costly event in a parallel environment, involving the coordinated control over a number of processors. This costliness, coupled with the relative rarity of inhibitor events, suggested a compromise solution: inhibitor events can be most easily handled by temporarily interrupting parallel processing altogether. The basic algorithm for our model was therefore adapted to follow a repetitive sequence of two basic steps: first, a random number of crystal association/dissociation steps is performed in parallel across all regions of the simulation space (the actual number of steps being based on the relative probabilities of crystal molecule and inhibitor molecule association/dissociation events); following this, parallel processing is interrupted to allow for a single inhibitor association/dissociation event. Experimental A series of whole crystal MC simulations were performed to investigate the inhibitory effects of AFP-I on hydrate sII formation. The unit cell of hydrate sII contains 136 water molecules, arranged so as to form 16 small and 8 large cages (Figure 1A). A simulation space structure was adopted that combined the spatial positioning and connectivity of the hydrate sII water molecules, with additional positions for gas molecules located at the center of each of the eight large cages, giving a total of 144 water and gas positions in the simulation space unit cell. Each of the eight central positions within the large cages was connected to the 28 water locations that comprise each large cage. The strength of each of these connections was set to 1/8th of the strength of inter-water location connections, reflecting the weaker dispersive forces acting between gas and water in hydrates than between waters. All simulations were run at 277K. As a first approximation, we adopted both the hydrogen bond energy value and the water association rate from our ice growth studies. In a first set of simulations that lacked inhibitor molecules, the gas association rate was varied to explore the effect of this parameter on hydrate growth rates and gas occupancy. The range of variance was between 0.625% and 100% of the water association rate. The model inhibitor, AFP-I, is a small alpha helical protein containing 38 amino acids (Figure  1B). For initial inhibitor studies, the simplifying assumption was made that the same AFP-I face that docks with ice Ih also docks with hydrate sII. This protein face is planar, with a repeating structure that should improve binding to crystalline structures (as is the case with ice), and its surface, comprised mainly of alanine residues, is hydrophobic, so that it reduces its solvent exposure upon binding. An initial set of inhibitors (ie. complex molecules) were built by first reducing the 136 water molecules in the hydrate sII unit cell to 7 distinct symmetric classes. One member from each of these classes acted as an anchor for AFP-I binding. AFP-I was positioned at each of the 12 vertices of a dodecahedron around the anchors, giving a three dimensional spectrum of docking orientations around the different classes of hydrate sII water molecules. 8 different docking models were created for each of these orientations by rotating the inhibitor 45º about an axis joining the anchor and the dodecahedron vertex, keeping the same AFP-I face oriented towards the anchor. This gave a total of 7 x 12 x 8 = 672 complex molecules for a first screening set, each representing a unique docking between AFP-I and a particular face of the hydrate sII. Volume sets were created for these inhibitors by determining which simulation space locations were carved out by each oriented docking. Neighbouring sets were created by discovering which simulation space locations were adjacent to the volume set locations. Uniform energetics values were assigned to these neighbouring locations. To standardize energetics between inhibitors, the energetics values assigned were scaled by the size of the energetics set (ie. by the number of neighbouring locations). A second set of inhibitors were built in an analogous manner, although this set contained only half as many inhibitors as only 4 different rotations about the anchor-vertex axis were used (0º, 45º, 90º, 135º). For this second set, different energetics were applied: neighbours adjacent to the binding face were given energetics values equivalent to ¼ of a hydrogen bond, whereas all other neighbours were given energetics values of 0. In effect, the energetics for the first set of inhibitors were uniform across inhibitor surfaces and therefore neutral to any particular AFP-I binding face, whereas those for the second set of inhibitors were given preferences for the selected binding face.  Gas Hydrate Simulations Without Inhibitors The accuracy of simulation results is entirely dependent upon the use of realistic association and  No. of Processors 1 2 4 8 16 24 48  Incremental Gainsb 1.00 0.92 0.81 0.65 0.47 0.34 0.17  Run Time 2668 1443 824 510 357 326 327  a  Identical simulation conditions (initial crystal size of 1 million molecules, run at 277K for 500 million simulation steps) were used for each run. b  Incremental gains measured by: Tserial / ( Tparallel * no. processors )  14  C ry s ta l S iz e  RESULTS AND DISCUSSION We have modified our existing MC crystal growth simulation software for studying ga hydrate growth in the presence of inhibitor molecules. Most significantly, the simulation software was parallelized to enable a screening for inhibitor docking models which exhibit significant hydrate growth inhibition. Parallelization was achieved by spatially partitioning the simulation space into symmetric regions and assigning one processor to manage crystal growth simulations in each region. Table 1 contains timing data from a series of runs of a test simulation, each using a different number of processors. All tests began with a 1 million molecule gas hydrate crystal, and growth proceeded for 500 million simulation steps in the presence of inhibitor. Performance improves substantially as the number of processors is scaled up to 16, although the incremental performance gain declines. Beyond 16 processors, performance gains are modest at best. Such performance degradation as more processors are available for parallel processing is not unexpected, since the incremental gains from added resources soon become offset by increased processing requirements for resource management. Specifically with regards to this application, a larger number of simulation space regions increases the percentage of locations that lie on regional boundaries, causing an increase in the amount of inter-process communication required for data structure management. All further simulations were done using 16 processors.  Table 1. Simulation run timesa by number of processors.  (M illio n s m o le c u le s )  Both sets of inhibitors were subjected to MC simulations. All simulations started with initial gas hydrate crystals containing approximately 500 000 water and gas molecules, and were run for 1 billion simulation steps. All software was written in the 'C' programming language. Simulations were run on the SunFire Cluster (a symmetric multiprocessor environment based on Sun's UltraSPARC processor, running the Solaris Operating Environment) at the High Performance Computing Virtual Laboratory (HPCVL) based in Kingston, Ontario. Parallelization was achieved using OpenMP.  12 0.63%  10  6.25% 12.50%  8  25%  6  50% 75%  4  100%  2 0 0  500  1000  1500  2000  2500  3000  S im u la tio n S te p (1 0 0 's m illio n s )  Figure 3. Gas hydrate growth as a function of the gas association rate (expressed as a percentage of the water association rate). dissociation rate parameters for all simulation molecules. As a first approximation, we adopted the water association and dissociation rate values from our previous studies of ice Ih crystal growth. These rates were determined by fine tuning the hydrogen bond energy and the association rate so as to achieve a series of stable ice crystals (crystals that neither grew nor melted over long simulation times) at a variety of temperatures, such that the radii of these crystals matched the critical radii for the simulation temperatures. Although the warmer simulation temperature used for gas hydrate growth simulations (277K versus 270-273K for ice growth simulations) will affect the hydrogen bond energy and water association rate somewhat, this effect is not anticipated to be large. With regards to the gas molecules, an initial series of simulations was  Table 2. Gas occupancy as a function of gas association rate. Gas Association Ratea  # Waters  # Gas  Ratio  Occupancy  6.25  3229  147  22.0  77.4%  12.5  591162  33444  17.7  96.2%  25  3205519  184854  17.3  98.0%  50  6290820  366414  17.2  99.0%  75  8033415  469503  17.1  99.4%  100  9092479  531920  17.1  99.5%  a  Gas association rate is expressed as a percentage of the water association rate.  performed to determine an appropriate range for the gas association rate parameter. A series of simulations was performed, all starting from a 100 000 molecule gas hydrate crystal and run for 3 billion simulation steps, with the gas association rate varied between 0.625% and 100% of the water association rate parameter. As seen in Figure 3, this had a dramatic effect on hydrate growth rates, providing us with a simple mechanism for tailoring growth as needed. Hydrate crystals were stable with the gas association rate set to 6.25% of the water association rate; lower values caused the hydrate to melt, while higher values caused it to grow. The water:gas ratio was also dependent on the gas association rate parameter (Table 2), ranging from 22.0:1 to 17.1:1 (amongst those crystals that did not melt) as the association rate increases. Although experimental results suggest that gas hydrate cage occupancy is near complete (full large cage occupancy yields a water:gas ratio of 17:1), we have elected to use somewhat lowered association rates to reduce the gas hydrate growth pressure in the preliminary growth inhibition simulations reported below. All crystals grown in these initial simulations had similar octahedral morphologies, an example of which is shown in Figure 4. Unlike in previous simulation work with ice Ih crystals where the crystal morphology could be altered by changes to the water rate parameters, we found the morphology of gas hydrate crystals to be highly resistant to change: our albeit limited exploration of parameter space always produce octahedral crystals, or no crystals at all.  Figure 4. Typical octahedral morphology of hydrate sII crystals from simulations. Water molecules are shown in red, and gas that are not completely entrapped are shown in gold. Gas Hydrate Simulations with Inhibitors Two screening procedures were performed to search for successful docking orientations between hydrate sII and AFP-I inhibitor. The first screening procedure consisted of 672 MC simulations, each involving a different docking orientation between gas hydrate and inhibitor (see Methods for details). Uniform energetics values scaled by the number of neighbouring locations (to give standardized aggregate energetics values across all inhibitors) were assigned to the connections between inhibitor and surrounding simulation space locations. As a result, no particular inhibitor face had preferred gas hydrate interactions over any other. The average number of simulation space locations occupied by these inhibitors, and the average number of neighbouring locations surrounding these inhibitors, was 178.7 and 295.8 respectively, with standard deviations of 5.2 and 20.5. The low standard deviation values for the statistics suggest that, while the orientations of these AFP-I inhibitors with respect to the underlying hydrate sII crystal vary dramatically, their spatial relationships to the crystal (ie. their volume and connectivity) vary little. All simulations in the first screening set used gas association rates set to 12.5% of the water association rate. Simulations began with 500 000 molecule gas hydrate crystals, and included a total of 1 billion simulation steps at a simulation temperature of 277K. We have performed a cursory examination  Figure 5. Various crystal morphologies from the first screening for gas hydrate inhibitors. Most crystals had octahedral morphology (A). A few had stunted growth at vertices (B). Others ranged from spherical (C) to cubic (D) to disk-like (E). Several showed oblong, uneven growth (F).  N u m b e r o f C ry s ta ls  250  (11.2)  200  (16.6) 150  (23.4) 100  (28.1)  (31.6) (35.5)  50  (36.2)  (0)  (39.3)  (44.2)  (51.1)  0%  0% 27  0-  29  0%  27 025  23  0-  25  0%  0% 21  0-  23  0% 19  0-  21  0%  19 17  0-  17  0% 015  13  0-  15  %  13 0-  90  11  -1  10  0% -9  0%  0  70  of the morphologies of the resulting 672 hydrate crystals. The great majority of these crystals had similar, roughly octahedral morphologies (Figure 5A). However, a spectrum of other morphologies was observed (Figure 5B-5F), ranging from octahedral with trimmed peaks, through spherical, cubic and disk-like. Several showed growth inconsistencies on one facet of their octahedrons, producing oblong non-symmetric shapes (Figure 5F). As expected, final inhibitor binding varied from crystal to crystal, both in location and in quantity. Figure 6 groups the final crystal sizes from these simulations, broken into ranges, along with the average number of bound inhibitors for crystals within each range. Crystal sizes are expressed in relation to the size of a reference gas hydrate crystal grown in the same simulation conditions without inhibitors. Unexpectedly, almost all of the inhibitor simulations produced hydrate crystals that were larger than the reference crystal, and some dramatically so. Moreover, there is a clear relationship between the average number of bound inhibitors and the final hydrate size. Rather than inhibit growth, these inhibitors appear to be enhancing it. A detailed examination of several of the simulation runs shows that the inhibitor molecules in these simulations remain bound to the gas hydrates for relatively short simulation times, despite having high numbers of hydrate contacts. We speculate that these inhibitors are functioning as gas hydrate nucleators: by forming transient interactions with the hydrate surface, these molecules are able to assist in the nucleation of new layers of hydrate growth at positions immediately adjacent to their binding sites, a feat attributable to their uniform energetics that do not actively dissuade adjacent growth (Figure 7). After new layers have been nucleated and are able to continue growth on their own, the inhibitors dissociate, leaving the newly nucleated layer free to fill in their vacated positions. A second screening was performed using a different set 336 inhibitors with altered energetics to determine whether stronger, longer lasting interactions between inhibitor and hydrate would promote growth inhibition. The binding strength between gas hydrate positions immediately adjacent to the AFP-I binding face were given strong, positive interaction values, while all other positions were  Ra n g e o f C r y s ta l S iz e s (e x p r e s s e d a s a % o f r e fe r e n c e c r y s ta l s iz e )  Figure 6. Number of crystals from the first screening, grouped by crystal size expressed as a percentage of the size of a reference crystal grown without inhibitors. The average number of bound inhibitors on the crystals in each range is given in parenthesis above each range. given interaction values of zero. In addition, we allowed inhibitors with a high degree of crystal contacts to achieve irreversible binding by setting the dissociation probability to zero once the binding strength achieved a minimal threshold (14 well placed connections to the hydrate crystal running along the length of the inhibitor binding face of AFP-I). Moreover, the gas association rate was lowered to 9.3% of the water association rate to further slow gas hydrate growth. Figure 8 shows the sizes of the final crystals from these simulations,  C ry s ta l S iz e  Figure 7. Schematic diagram of proposed hydrate nucleation by inhibitors in simulations. Top left: an inhibitor first binds to the clathrate surface. Middle: positive energetics between the inhibitor faces adjacent to the binding face and the hydrate nucleate a new hydrate layer. Bottom right: the inhibitor dissociates, leaving the nucleated layer free to grow.  N u m b e r o f C ry s ta ls  100  (M illio n s m o le c u le s )  2  1 .6 65-75% 75-85%  1 .2  85-95% 125-135%  0 .8  165-175% R eferenc e  0 .4  0 0  200  400  600  800  1000  S im u la tio n S te p s (1 0 's m illio n s )  Figure 9. Growth rates of representative hydrate crystals from the second inhibitor screen. Crystals were chosen from several of the different crystal size ranges shown in Figure 8. Growth inhibition was almost achieved in the best simulation (pink).  (120) (237)  (132)  75  (282) 50  (153) 25  (298)  (140)  (113) (104)  (0)  (40)  5%  5% 16  5-  17  5%  16 515  14  5-  15  5%  5% 13  5-  14  5%  13  12  5-  511  12  %  5% 10  5-  11  05  5%  -1  -9  95  -8  5% 85  75  65  -7  5%  0  Ra n g e o f C r y s ta l S iz e s (e x p r e s s e d a s a % o f r e fe r e n c e c r y s ta l s iz e )  Figure 8. Number of crystals from the second screening, grouped by crystal size expressed as a percentage of the size of a reference crystal grown without inhibitors. The average number of bound inhibitors on the crystals in each range is given in parenthesis above each range. again expressed in relation to a reference crystal grown in the same simulation conditions without inhibitors. Substantially less of the resulting crystals were larger than the reference crystal in this second screening as compared to the first, and most of those that were larger were only modestly so. Moreover, there were a number of simulations that actually reduced gas hydrate growth. Figure 9 shows plots of hydrate growth over simulation time for a number of  Figure 10. Various crystal morphologies from the second screening for gas hydrate inhibitors. Most crystals had much higher inhibitor coverage, ranging from a majority of the surface (A) to full coverage (C). Octahedral morphology (A) was dominant, though some crystals had spherical (B) and cubic (C) morphologies. Numerous had non-symmetric shapes (D), reflecting uneven hydrate growth in response to inhibitor binding. Many also had plane-specific inhibitor binding (E), which was unexpected because of the inherent symmetries in the hydrate sII unit cell. (F) depicts vertexspecific inhibitor binding, and the resulting growth inhibition at these locations.  representative crystals from various growth ranges. Clearly, the increased binding strength of the inhibitor binding face, together with reduced binding strengths of all other inhibitor faces, has had dramatic effects on gas hydrate growth. Interestingly, although the volumes and energetics of the inhibitors were similar, they produced a wide spectrum of effects upon hydrate growth, with resulting crystals ranging between 0.7 and 1.66 times the size of the reference crystal. This wide variation suggests that the specific docking arrangement is also critical to inhibitor activity. Similar conclusions were drawn from work with ice crystal growth in the presence of inhibitors [14, 15]. From a cursory inspection of the resulting crystals, it is clear that far more inhibitor remains bound to the crystal surface than was observed in the first inhibitor screens. Figure 10 shows a representative final crystal. CONCLUSIONS We have adapted our novel Monte Carlo whole crystal simulation technique for studying gas hydrate sII crystal growth, both in the presence and absence of model inhibitors based on AFP-I. We have performed a number of inhibitor screening studies involving hundreds of potential hydrate – AFP-I docking orientations in an attempt to identify inhibitor and docking characteristics that promote gas hydrate growth inhibition. Our preliminary inhibitor screenings suggest that there is a connection between binding energetics and gas hydrate formation rates. Our first inhibitor screen unexpectedly resulted in an almost universal increase in hydrate growth for all inhibitors, and an examination of several of the simulations revealed two possible reasons for this phenomenon, both related to the inhibitor energetics. First, the uniform nature of the energetics of these inhibitors meant that neighbouring locations on the sides of the inhibitor directly adjacent to the binding plane would gain some additional stability from the inhibitor. This extra stability was enough to tip the association/dissociation balance in favour of association, effectively promoting the nucleation of new hydrate growth planes on the hydrate surface. Secondly, the inability of the bound inhibitors to form long lasting connections with the hydrate crystals seriously compromised their abilities to inhibit growth. The combination of these two  properties effectively transformed the inhibitors in this first screen into nucleators: inhibitor binding resulted in increased nucleation at immediately adjacent locations; subsequently, inhibitor dissociation enabled newly nucleated growth planes on the gas hydrate surface to complete growth without impingement. The importance of these energetics features was confirmed with a second screening process using a set of inhibitors with energetics that reflected (i) stronger interactions to hydrate positions connected to their binding faces, (ii) neutral interactions to all other neighbouring positions, and (iii) irreversible binding when sufficient binding surface complementarity is achieved. While the majority of simulations involving inhibitors in this second screening still failed to display substantial gas hydrate growth inhibition, many did reduce growth, with several producing substantially smaller crystals. The results of both screening tests – nucleation in the first set and inhibition in the second set – showed significant deviation between individual inhibitors. We conclude from this that, in addition to the identified energetics features listed above, both nucleation and inhibition are dependent upon specific binding orientations. This is consistent with previous findings on ice growth inhibition [14, 15]. We close by stressing that these simulations are preliminary, and much work remains to be done. Indeed, this work stimulates more questions than it answers: is irreversible binding absolutely necessary for inhibition, or can inhibition be achieved while maintaining a basal level of inhibitor dissociation? Which specific features of the energetics promotes gas hydrate nucleation? Which ones promote inhibition? What are the characteristics of the binding orientation best able to increase or suppress hydrate growth? We are turning to these and other specific questions in the next phase of our research. We believe that answers to these questions will lead to a better understanding of the mechanisms of gas hydrate growth and growth inhibition. ACKNOWLEDGEMENTS This research is supported by the High Performance Computing Virtual Laboratory in Kingston, Canada. BW is a recipient of an NSERC graduate studentship, and the Sun Microsystems of Canada Scholarship in Computational Sciences and  Engineering. VKW is a recipient of an NSERC grant. ZJ is a Canada Research Chair in Structural Biology. We are thankful to Jo Willan for her assistance in crystal screening. REFERENCES [1] Kvamme B, Kuznetsova T, Aasoldsen K. Molecular dynamics simulations for selection of kinetic hydrate inhibitors. J. Mol. Graphics Modelling 2005; 23:524-536. [2] Sloan ED Jr. Hydrate Engineering Monograph Vol 21. Richardson, Texas: Society of Petroleum Engineers, 2000. [3] Gao S, House W, Chapman WG. Detecting Gas Hydrate Behavior in Crude Oil Using NMR. J. Phys. Chem. B 2006; 110:6549-6552. [4] Sloan ED Jr. Clathrate Hydrates of Natural Gases. 3rd ed.; New York: Crc Press Llc., 2007. [5] Schicks JM, Naumann R, Erzinger J, Hester KC, Hoh CA, Sloan ED Jr. Phase Transitions in Mixed Gas Hydrates: Experimental Observations versus Calculated Data. J. Phys. Chem. B 2006; 110:11468-11474. [6] Patchkovskii S, Tsu JS. Thermodynamic stability of hydrogen clathrates. Proc. Natl. Acad. Sci. USA 2003; 100:14645-14650. [7] Konstantin AU, Lu H, Enright GD, Ratcliffe CI, Ripmeester JA, Chapman NR, Riedel M, Spence G. Single Crystals of Naturally Occurring Gas Hydrates: The Structures of Methane and Mixed Hydrocarbon Hydrates. Angew. Chem. Int. Ed. 2007; 46:8220-8222. [8] Moon C, Hawtin RW, Rodger PM. Nucleation and control of clathrate hydrates: insights from simulation. Faraday Discuss. 2007; 136:367-382. [9] Zeng H, Wilson LD, Walker VK, Ripmeester JA. Effect of Antifreeze Proteins on the Nucleation,  Growth, and the Memory Effect during Tetrahydrofuran Clathrate Hydrate Formation. J. Am. Chem. Soc. 2006; 128:2844-2850.  [10] Moon C, Taylor PC, Rodger PM. Molecular Dynamics Study of Gas Hydrate Formation. J. Am. Chem. Soc. 2003; 125:47064707. [11] Alavi S, Ripmeester JA, Klug DD. Molecular-dynamic Simulations of Binary Structure II Hydrogen and Tetrahydrofuran Clathrates. J Chem. Phys. 2006; 124:1470414709. [12] Alavi S, Ripmeester JA, Klug DD. Molecular Dynamics Study of the Stability of Methane Structure H Clathrate Hydrates. J. Chem. Phys. 2007; 126:124708-124713. [13] Anderson BJ, Tester JW, Borghi GP, Trout BL. Properties of Inhibitors of Methane Hydrate Formation via Molecular Dynamics Simulations. J. Am. Chem. Soc. 2005; 127:17852-17862. [14] Wathen B, Kuiper M, Walker VK, Jia Z. A New Model for Simulating 3-D Crystal Growth and its Application to the Study of Antifreeze Proteins. J. Am. Chem. Soc. 2003; 125:729-737. [15] Wathen B, Kuiper M, Walker VK, Jia Z. New Simulation Model of Multicomponent Crystal Growth and Inhibition. Chem. Eur. J. 2004; 10:1598-1605. [16] Gilmer G. Computer Models of Crystal Growth. Science 1980; 208:355-363.  


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items