Solid-liquid Equilibrium and Free Energy

1
Solid-liquid Equilibrium and Free Energy
Calculation of Hard-sphere, Model protein, and
Lysozyme Crystals
Jaeeon Chang
Center for Molecular and
Engineering Thermodynamics Department of
Chemical Engineering University of Delaware
2
Overview

• Determination of liquid-solid transition using
  histogram reweighting method and expanded
  ensemble MC simulations
• Fluid-solid phase equilibria of patch-antipatch
  protein model
• The combined simulation approach of atomistic
  and continuum models for the thermodynamics of
  lysozyme crystals

3
Objective
To develop a generic method to predict
liquid-solid phase equilibria from Monte Carlo
simulations. Motivation Methods involving
particle insertion scheme used in Gibbs ensemble
MC are not applicable to dense liquids and
solids. The previous method based on equation of
state requires fitting of simulation data to an
assumed form of EOS. Kofkes Gibbs-Duhem
integration method requires one known point
on the coexistence curve.
4
Hard sphere fluid and fcc crystal (Test model)

• The simplest nontrivial potential model
• no vapor-liquid transition, athermal
  solid-liquid transition
• Reference system for perturbation theory
• Model for colloid system

• Canahan-Starling equation of state for liquid
  phase

5
Outline of the Methodology
? Basic principle Tliq Tsol pliq
psol ?liq ?sol
? A series of NPT MC simulations are performed
to separately construct equations of state for
liquids and solids using histogrm reweighting
method
? To obtain the chemical potentials of liquid
branch an accurate estimate at a particular
density should be provided either from direct
simulations (Widom method, Free energy
perturbation method, Bennett acceptance method)
or from the integration of the equation of state
from from zero density to the liquid density.
? For the chemical potentials of solid branch,
the free energy at a particular density is
obtained using Einstein crystal and the expanded
ensemble method.
6
Histograms from NPT Monte Carlo simulations
? Probability for a single histogram reweighting
h
? Hard spheres
V, E
7
Histogram Reweighting Method in NPT ensemble
? Composite probability for multiple histograms
(not normalized)
Cis are determined in a self-consistent
manner
? Average properties
? Chemical potential
A known value of free energy is required to
specify C0
8
Construction of Equation of State for Hard spheres
9
Free Energy of Solid
The classical Einstein crystal as a
reference Variation of potential from the
reference to the system of interest
Einstein crystal (reference)
Repulsive core turned on
Einstein field turned off
The expanded ensemble method hopping over
subensembles
Ref) Lyubartsev et al., J. Chem. Phys. 96, 1776
(1992). Chang and Sandler, J. Chem. Phys.
118, 8390 (2003).
10
Free Energy of Hard-sphere Fcc Solid
11
Determination of Solid-liquid Transition

• Equilibrium properties
• p 11.79, ?liq 0.944, ?sol 1.045, ?
  16.30
• c.f) Lee and Hoover (1968)
• p 11.70, ?liq 0.943, ?sol 1.041

12
Conclusions

• Using the combination of the histogram
  reweighting and expanded ensemble simulation
  methods a new generic algorithm for predicting
  liquid-solid equilibria is proposed.
• The liquid-solid equilibria for hard-sphere
  systems of varying size up to 1372 particles are
  studied, and the limit for the infinitely large
  system is accurately determined.

13
Overview

• Determination of liquid-solid transition using
  histogram reweighting method and expanded
  ensemble simulations
• Fluid-solid phase equilibria of patch-antipatch
  protein model
• The combined simulation approach of atomistic
  and continuum models for the thermodynamics of
  lysozyme crystals

14
Introduction

• Background
• 3-D protein structure by X-ray crystallography
• Crystallization windows correlated with slightly
  negative values of B2
• L-L equilibria described by isotropic
  short-range interaction
• Anisotropic model necessary for F-S equilibria

• Objectives of this work
• Understanding of the role of anisotropic
  interactions on fluid-solid equilibria of protein
  solutions from exact computer simulations
• Comparison of the free energies of different
  crystal structures and phase diagram involving
  multiple solid phases

15
Patch-antipatch potential model of globular
proteins
Three patch-antipatch pairs in perpendicular
directions A narrow range of orientation
ap 12º
45 for chymotrypsinogen
Ref Hloucha et al., J. Crystal Growth, 232, 195
(2001)
16
Crystal structures
High T Low T Low T
simple cubic SC
orientationally disordered face centered
cubic FCC(d)
orientationally ordered face centered
cubic FCC(o)
17
Free energy (Fluid)
Thermodynamic integration method over the
equation of state
18
Free energy of the Einstein crystal
Constraining potential
Free energy
19
Free energy (Solid)
The classical Einstein crystal as a
reference Variation of potential from the
reference to the system of interest
Einstein crystal (reference)
Repulsive core turned on
Attractions turned on
Einstein field turned off
The expanded ensemble method a direct measure
of free energy difference
Ref) J. Chang and S. I. Sandler, J. Chem. Phys.
118, 8390 (2003)
20
Coexistence at a fixed temperature
T 1.0, ?p 1 (isotropic model)
21
Coexistence for isotropic model
Density
Osmotic pressure
22
Coexistence for anisotropic model with ?p 5
Density
Osmotic pressure
The ordered phases SC and FCC(o) are stable
at low temperatures
23
Phase diagram for anisotropic model with ?p 5
Chemical potential of solid beyond saturation
pressure
Saturated solids
Phase diagram
FCC(o) cannot be prepared directly from solutions
24
Conclusions

• More than one crystal structures are compatible
  even with the simple anisotropic protein model.
• Phase diagram involving multiple solid phases are
  determined from the histogram reweighting method.

25
Overview

• Determination of liquid-solid transition using
  histogram reweighting method and expanded
  ensemble simulations
• Fluid-solid phase equilibria of patch-antipatch
  protein model
• The combined simulation approach of atomistic
  and continuum models for the thermodynamics of
  lysozyme crystals

J. Phys. Chem. B, 109, 19507 (2005)
26
Introduction

• Protein crystal
• Protein crystals are used to determine
  three-dimensional structure of proteins.
• Trial-and-error screening methods to find
  crystallization conditions
• Weak attractive interactions for the
  crystallization to occur

Hen egg white lysozyme (193L)

• Lysozyme
• Natural antibiotic enzyme to break the cell wall
  of bacteria
• Several crystal forms depending upon the solution
  conditions ( T, pH, Ionic strength and species)
• Solubility, heat of crystallization are known
  experimentally.
• Phase transition between tetragonal (low T) and
  orthorhombic (high T) forms occurs near the room
  temperature

27
Objective and approach

• To compare the thermodynamic properties of
  tetragonal and orthorhombic crystals of hen egg
  white lysozyme

Crystallographic structure
Boundary element method
NVT Monte Carlo simulation
Elec. A
Expanded ensemble MC simulation
Gibbs-Helmholtz relation
vdw U
Elec. A
Elec. U, S
vdw A
Thermodynamic properties
28
NVT Monte Carlo simulations

• Potential model (implicit water)
• Semi-empirical model of Asthagiri et al.,
  Biophys. Chem. (1999)
• 0.5 ? OPLS force field, r lt 6Å
• Hamaker interactions ( H 3.1 kT), r gt 6Å
• System
• 16 protein molecules (rigid body)
• NVT MC at experimental density and at 298 K
• Translation 0.1Å, Rotation 1º

orthorhombic form (PDB 1F0W)
tetragonal form (PDB 193L)
29
Free energy of crystal from expanded ensemble MC
Variation of potential from the reference to the
system of interest
Einstein crystal (known free energy)
Einstein field turned off Full interaction
? 0 0 lt ? lt 1 ? 1
The expanded ensemble MC method hopping over
subensembles
Ref) Lyubartsev et al., J. Chem. Phys. 96, 1776
(1992). Chang and Sandler, J. Chem. Phys.
118, 8390 (2003).
30
Boundary element method for electric potential

• Protein domain dielectric const. ? 4

• Aqueous domain ( a single protein in solution )
  dielectric const. ? 80

• Aqueous domain of the Crystal is in the Donnan
  equilibrium with the solution.
• Optimally linearized Poisson-Boltzmann
  equation

• Electrostatic free energy

31
Boundary elements and charge distribution in the
protein

• 1,000-2,000 Triangular elements were obtained
  using Connollys program with subsequent
  simplification procedure.
• The charges are placed at 30 ionizable residues
  (Asp, Glu, Arg, His, Lys, Tyr) and the C and N
  termini.
• The Henderson-Hasselbalch equation using
  experimental pKa data.

32
VDW energy from MC simulations

• The VDW interactions in the tetragonal form are
  less attractive than in the orthorhombic form.
• Noticeable variations in the energy among the PDB
  structures are observed due to the variations in
  the side chain conformations.
• The Lennard-Jones interactions within 6 Å
  dominate over the water-mediated Hamaker
  interactions at longer distances.

33
Free energy of lysozyme crystal (pH 4.5 and I
0.36 M)

• Standard state at 1mol/L A/NkT ln(?3/1660)
  1
• The experimental values are close to each other
  since the transition occurs near 298 K.
• Electrostatic contribution to the free energy is
  repulsive.
• The predicted Helmholtz energies are less than
  the experimental values.

Tetragonal
Orthorhombic
34
Energy of crystallization (pH 4.5 and I 0.36 M)

• The VDW contributions play a dominant role in
  both crystals.
• For the tetragonal crystal, the predicted energy
  is acceptable considering the wide range of the
  reported experimental values from 40 to 140
  kJ/mol.
• The less attractive energy from experiment
  suggests energetically unfavorable release of
  water molecules from crystal contacts.
• The larger disagreement for the orthorhombic
  crystal form indicates a large difference in the
  solvation structure.

Tetragonal
Orthorhombic
water release
35
Entropy of lysozyme crystal (pH 4.5 and I 0.36
M)

• The VDW entropy for the tetragonal form is in
  good agreement with the mean field theory.
• The electrostatic contributions to the entropy
  are negative, arising from the reorganization of
  water molecules and ions.
• For the tetragonal form, there should be a
  release of about 4 water molecules upon
  crystallization ( the entropy change on
  the melting of ice 22 J/mol/K).

water release
Tetragonal
Orthorhombic
36
Experimental evidence of distinct hydration
structures

• A water molecule is counted as a bridging
  molecule if it is also close to another protein.
• Whereas the total number of hydrated water
  molecules is almost the same, there is a decrease
  of about ten bridging water molecules for each
  protein in the orthorhombic crystal forms.
• Additional water molecules are expelled from
  between the contacting surfaces when a lysozyme
  molecule becomes part of an orthorhombic crystal,
  which is an energetically less favorable but
  entropically more favorable process.

gt
37
Conclusions

• We have carried out Monte Carlo simulations of
  the hen egg white lysozyme crystals at the
  atomistic level and the boundary element
  calculations to solve the Poisson-Boltzmann
  equation for the electrostatic interactions.
• The crystallization energy of the tetragonal
  structure agrees reasonably well with
  experimental data, while there is a
  considerable disagreement for the orthorhombic
  form.
• A large difference in the experimental energy of
  crystallization between the two crystals
  indicates energetically unfavorable solvation in
  the orthorhombic form.
• The much higher value of the entropy of the
  orthorhombic crystal is explained in terms of the
  entropy gain of the water molecules released
  during the crystallization.

38
Acknowledgements

• National Science Foundation
• Department of Energy
• Prof. Stanley Sandler
• Prof. Abraham Lenhoff
• Dr. Jeffrey Klauda (NIH)
• Dr. Stephen Garrison (NIST)
• Mr. Gaurav Arora