Transport and Structure in Nanoscale Channels by Gregory W i l l i a m Lakatos A THESIS S U B M I T T E D IN PARTIAL F U L F I L M E N T O F THE REQUIREMENTS FOR T H E DEGREE OF DOCTOR OF PHILOSOPHY in T H E F A C U L T Y OF G R A D U A T E STUDIES (Physics) T H E UNIVERSITY OF BRITISH C O L U M B I A June 21, 2006 © Gregory W i l l i a m Lakatos, 2006 11 Abstract Driven by the rapidly advancing fields of nano- and biotechnology, there has been an explosion of interest i n molecular transport and structure formation on small length scales. A canonical model for the transport of particles along one dimensional pathways i n nanoscale channels is the Totally A s y m m e t r i c Simple Exclusion Process ( T A S E P ) . After introducing the standard T A S E P , modifications of the T A S E P designed to increase its utility in modeling biological transport processes are described. One variant of the T A S E P is particularly suitable for modeling protein translation, and the results of using this variant to investigate the effects of slow-codons on the translation process are discussed. A related topic is the voltage-driven translocation of D N A hairpins through membraneembedded nanopores. M o t i v a t e d by recent experiments, a stochastic model is developed that couples the translocation and dehybridization of the D N A hairpin. T h i s model is used to explore the behaviour of the mean translocation time of hairpins as a function of driving voltage, and two translocation mechanisms are identified and discussed. Finally, the adsorption and equilibrium structures of water i n the interior of ion-bearing nanoscale pores are considered. T h e behaviour of water and ions under confinement is critical to the functioning of biological ion channels and nanoporous filters. Here, the adsorption isotherms of water are examined, and the layered structures formed by the confined water are described. iii Contents Abstract ii Contents iii List of Tables v List of Figures vi Preface xiv Acknowledgements xvii 1 Introduction 1 2 T h e Totally A s y m m e t r i c Simple Exclusion Process 4 3 Totally A s y m m e t r i c Exclusion Processes w i t h F i n i t e Size Particles 3.1 3.2 3.3 3.4 4 Background M e a n F i e l d Theories 3.2.1 Simple M e a n F i e l d Methods 3.2.2 Refined M e a n F i e l d Theory for the M a x i m a l Current Phase 3.2.3 Refined M e a n F i e l d Theory for Boundary-limited Currents 3.2.4 M a t c h i n g M e a n F i e l d Phase Boundaries Monte Carlo Simulations 3.3.1 Currents 3.3.2 Particle Densities 3.3.3 Phase Boundaries Summary . . . . . . . A Totally A s y m m e t r i c Exclusion Process with Periodic Structure . . . 4.1 Background 4.2 M o d e l and Methods 4.3 M e a n F i e l d Theories 4.3.1 Simple M e a n F i e l d Methods 4.3.2 Refined M e a n F i e l d Methods 4.4 Monte Carlo Simulations 4.4.1 Currents 4.4.2 Densities 4.4.3 Phase D i a g r a m . 12 12 13 13 17 19 22 23 24 27 29 33 34 34 35 36 36 38 48 48 52 56 Contents 4.5 4.6 iv Other Periodic Arrangements Summary 56 61 5 Clustered Bottlenecks in m R N A Translation and P r o t e i n Synthesis . . 5.1 Background 5.2 M o d e l 5.3 Results and Discussion 5.4 Summary 63 63 65 69 74 6 First Passage T i m e s of D r i v e n D N A H a i r p i n U n z i p p i n g 6.1 Background 6.2 M o d e l 6.3 Results and Discussion 6.4 Summary 76 76 78 86 95 7 8 Water A d s o r p t i o n in Ion bearing Nanopores 7.1 Background 7.2 M o d e l and Simulation Methodology 7.2.1 Molecular Models 7.2.2 Simulation Methodology 7.3 Results and Discussion 7.3.1 A d s o r p t i o n Isotherms 7.3.2 F i l l e d State Structures 7.3.3 A d s o r p t i o n In Ion-Bearing Pores 7.3.4 Pressures 7.3.5 G r a n d Potentials 7.3.6 Boundary Effects 7.4 Summary • • 97 97 98 98 104 110 110 114 116 132 134 139 143 : Conclusions and Suggestions for Future W o r k 146 Bibliography A A n implementation of the rithm 152 Bortz-Kalos-Lebowitz Monte Carlo Algo158 List of Tables 6.1 6.2 6.3 6.4 7.1 T h e fixed size and rate parameters used i n the two dimensional stochastic model. T h e size parameters are set by the structure of the a-hemolysin pore and the D N A strand, and are not used as fitting parameters List of fitting parameters used i n the two dimensional stochastic model. . D N A hairpin sequences [52] along w i t h the hybridization energies predicted by M F O L D ( T = 288) and the 2-3 free energy model. For more detail see [4]. Bases i n the hairpin hybridization region are i n bold, the defect bases are underlined, and the other bases form the loop segment of the hairpin. . Parameters from fitting to the experimental results for the mean escape time of D N A hairpins through an a-hemolysin pore of M a t h e et al. [52]. T h e fit parameters sets were used to generate the mean escape time voltage curves shown i n F i g . 6.3(a) T h e potential parameters (at T = 298K) used for the different species [60, 96, 104]. A l l Lennard-Jones parameters are for interactions between the listed species and S P C / E water. Cross-interaction parameters were computed using Lorentz-Berthelot mixing rules [105]. T h e carbon parameters were used for the wall potential given i n E q . 7.2 79 85 87 88 99 vi List of Figures 2.1 3.1 3.2 3.3 3.4 (a) T h e standard totally asymmetric simple exclusion process ( T A S E P ) . Particles are injected on the left at a rate a, traverse the lattice from left to right by making hops with an attempt rate p, and are removed from the lattice w i t h an attempt rate (3. N o more than one particle can occupy a lattice site at a time, (b) T h e three solution regimes of the standard T A S E P . Note that here p has been set equal to 1. I N S E T S : T h e density profile along the T A S E P lattice. In the entry and exit limited phases two density profiles are possible depending on the values of (a,/3) relative to the a = 1 — P line (the dashed line i n the plot) E n t r y and exit mechanisms, (a) P a r t i a l entry (b) complete entry into the first d sites (c) partial exit (d) complete exit as the particle reaches the last site (a) A typical configuration of particles associated w i t h the calculation of Z(n, L) given by E q . 3.11. The left edge of the particles is noted by a black dot, while the exclusion zone associated w i t h each particle is shown with a blue rectangle (b) A typical configuration for the calculation of Z(n, L — 1). State enumerations for particles of size d = 4 near the lattice boundaries, (a) T h e d+l distinct states at the entrance region, (b) T h e d + 1 distinct states of the exit region. Note that a 0 indicates the absence of the left edge of a particle at the lattice site. Thus i n the states displayed in (b), all the 0 sites may be occluded by part of a particle, but cannot carry a left edge. In b o t h (a) and (b), all the states that are w i t h i n d moves of state Si are included T h e currents w i t h i n each of the solution regimes as a function of particle size (d). Results from simple mean field, refined mean field, and Monte Carlo simulations are shown, (a) T h e entry limited current JL as a function of d for (3 = 10 and a = 0.1. (b) T h e m a x i m a l current J for a = ft = 10. (c). T h e exit limited current JR for a = 10 and (3 = 0.1 Currents as a function of a or j3 for d = 2 ,4, and 8. (a) T h e entry limited to m a x i m a l current profile w i t h (3 = 10. (b) T h e exit limited to maximal current profile w i t h a = 10. T h e points are the M o n t e Carlo results, the solid lines are the refined mean field predictions ( E q . 3.24), and the dashed lines are the simple mean field predictions (Eq. 3.8) Boundary densities i n each of the three current regimes. T h e solid lines show the predicted densities from the refined mean field approximations (Eq. 3.24) 6 15 17 19 m a x 3.5 3.6 25 26 28 List of Figures 3.7 3.8 3.9 4.1 vii Representative scaled density (aid) profiles i n the three current regimes, (a) T h e entry rate limited regime (a = 0.1, (3 = 1). (b) T h e m a x i m a l current regime (a = (3 = 1). (c) T h e exit rate limited regime (a = l,/3 = 0.1). T h e boundary regions are shown i n the insets. T h e strong splitting of the densities into high density and low density branches, is a product of the close-packing of the particles i n the exit limited regime, and the fact that particle positions are determined by the location of the particles' left edges. 30 T h e phase diagram for the T A S E P w i t h d = 3 particles derived from the refined mean field results (Eq. 3.24) and the simple mean field results (Eq. 3.8). T h e points are the transitions estimated from M o n t e Carlo simulations. 31 (a) T h e discontinuity i n the second derivative of the steady-state current w i t h respect to a or B at the transition between the m a x i m a l current and the entry limited regimes. T h e simple mean field approach (Eq. 3.9) predicts an increasing discontinuity w i t h d, while the refined mean field approach ( E q . 3.25) predicts a j u m p that decreases w i t h d. (b) T h e critical values a* and B* at which the boundary limited to m a x i m a l current transitions occur (Eqs. 3.8 and 3.24) 32 T h e periodic dual-rate totally asymmetric exclusion model. A t every T lattice sites, the particle movement rate is p . A t all intervening sites, the movement rate is p\. A l l other aspects of the model are identical to those of the standard T A S E P T h e arrangement of lattice sites used for an L = 4 finite segment mean field ( F S M F ) approximation to the T = 2 T A S E P . T h e master equation for the state of the four marked lattice sites is solved exactly. T h e four sites are coupled to the rest of the lattice by assuming an effective injection rate of « e f f = P2&4 on the left, and an effective extraction rate of B $ = ^2(1 — en) 2 4.2 35 e 4.3 on the right T h e steps i n the algorithm to generate the transition m a t r i x for a T A S E P . For the purposes of illustration a three site lattice has been used. Each possible occupancy of the lattice is treated as a bit pattern, and each state is labelled w i t h the corresponding decimal value (i.e. lattice occupancy O i l —» state 3). T h e states are divided into two groups; states where the first lattice site is occupied (1-states), and states where the first lattice site is empty (0-states). Regardless of the number of lattice sites i n the T A S E P , the transitions between the two classes of states always occur between the first half of the 1-states where the second lattice site is unoccupied and the second half of the 0-states where the second lattice site is occupied, (solid arrows i n the figure). Calling the algorithm recursively on b o t h the 1states and the 0-states, and ignoring the highest order bit (i.e. the leftmost lattice site), enumerates all the remaining state transitions (dashed arrows). F i n a l l y the transitions between each 0-state and each 1-state generated by injection at the left edge of the lattice (dotted arrows) are enumerated. . . 43 44 List of Figures 4.4 4.5 viii (a) T h e M o n t e Carlo ( M C ) , simple mean field, refined mean field, and L = 12 finite-segment mean field ( F S M F ) predictions for the T = 2 T A S E P current in the m a x i m a l current phase (a = 10, /? = 10). (b) M a x i m a l currents for dual-rate T A S E P s of various periods. T h e F S M F estimates shown by solid lines were produced using a single period of the lattice as the finite segment. W h i l e providing reasonable estimates, the deviation between the M C (a = 10, B = 10) and F S M F (lines) results increases as the period (T) increases. However (c) displays the percent difference (£) between the M o n t e Carlo and F S M F results for a T = 5 T A S E P , and shows that the increase i n the F S M F error w i t h increasing T can be partially reduced by increasing the segment length (a) Current profiles along the a direction i n the T — 2 dual-rate T A S E P phase plane (B — 10). T h e points are Monte Carlo results, the dashed lines were produced using E q . 4.4, and the solid lines were produced using an N = 8 F S M F approach. E q . 4.20 and the F S M F approximation are indistinguishable, (b) Current profiles for dual-rate T A S E P s of various periods w i t h p = 0.25. T h e F S M F results (lines) were produced using an L = T F S M F approach (/? = 10) M o n t e Carlo density profiles for a T = 2 T A S E P w i t h p = 0.1. Comparison w i t h F i g . 2.1 shows that the profiles in the T — 2 T A S E P retain the same qualitative character as the profiles i n the standard T A S E P . Note that in the boundary limited regimes the density profiles are approximately constant near the rate limiting boundary, consistent w i t h the assumptions leading to Eqs. 4.18 and 4.19 M a x i m a l current phase (a = 10,/? = 10) densities for a T = 2 T A S E P from Monte Carlo simulations, the refined m a x i m a l current mean field method (Eq. 4.14), the simple mean field method (Eq. 4.2), and an L — 4 finite segment mean field approach. T h e Monte Carlo results are averages over the central 50 sites of a 2000 site lattice. T h e simple mean field assumption produces the largest errors at small values of p where the density correlation between adjacent sites is expected to be large 49 2 4.6 4.7 51 2 52 2 4.8 53 T h e bulk density profiles for the T = 5 and T = 9 dual-rate T A S E P s w i t h p2 = 0.2, i n the m a x i m a l current (a = 10,0 = 10), entry rate limited (a = 0.05, 0 = 10), and exit-rate limited (a = 10,0 =. 0.05) phases. The plots show the good agreement between the M o n t e - C a r l o (points) and L = T F S M F (lines) results. T h e M C results are the central 5 periods (T = 5), and central 3 periods (T = 9) of the lattice 54 ix List of Figures 4.9 T h e density profile along the a direction i n the T = 2 dual-rate T A S E P phase plane ((3 = 10). T h e dashed lines were produced using E q . 4.4, while the solid lines were produced by an L = 8 finite-segment mean field approach. T h e M o n t e Carlo results are averages over the central 50 sites of the T A S E P lattice. T h e F S M F results are indistinguishable from those of E q . 4.20. A l t h o u g h the solution for cr is the same i n Eqs. 4.4 and 4.20, the predicted value of a* differs between the two mean field theories. T h i s is the reason for the significant difference between the o profiles predicted 2 2 by the two mean field methods at small values of p 4.10 (a) T h e phase diagram for the T = 2 two-rate T A S E P . T h e dashed lines are produced by the simple mean field theory ( E q . 4.7). T h e dotted lines are produced by the m a x i m a l current results (Eqs. 4.16, 4.17), and the solid lines are produced by the boundary-limited results(Eqs. 4.20, 4.21), while the points are the Monte Carlo results, (b) T h e phase diagram for the two-rate T A S E P at various periodicities, w i t h p = 0.1. To simulation accuracy, the transition between the high and low density phases occurred along the a = (3 line for all the periodicities displayed. Points are Monte Carlo results, lines are drawn as a guide to the eye 2 55 2 57 4.11 Variations on the periodic dual-rate T A S E P . T h e three lattices displayed correspond to (a) 8 — 1, (b) 8 = 2, (c) 8 = 3. In each case the lattice length was chosen to ensure the resulting T A S E P possessed particle-hole symmetry. 58 4.12 Phase diagrams for the T = 3 dual-rate T A S E P for various values of 8 and N derived from Monte Carlo simulations (p = 0.25). T h e lines are the best-fit phase boundaries for T A S E P s w i t h values of N that satisfy E q . 4.28. T h e points are the Monte Carlo predictions for the phase boundaries of T A S E P s w i t h values of N that violate E q . 4.28. Note that i n b o t h the (8 = 1,/V = 1203) T A S E P and the (8 = 2,N = 1204) T A S E P the last defect is located at lattice site N — 2 and the T A S E P s share a common value for /?*. Additionally, T A S E P s (8 = 2, N = 1203), (8 = 3,N = 1204), and (8 = 1, N = 1205) all share a common location for the last defect site, 2 and share a common value of (3* 5.1 60 T h e m R N A translation/protein synthesis process. Ribosomes move unidirectionally along the m R N A as t R N A s (not shown) deliver the appropriate amino acid to the growing protein. Codons w i t h low concentrations of corresponding t R N A ("slow-codons") result i n bottlenecks that locally slow ribosome motion. A totally asymmetric simple exclusion process on N lattice sites is used to model the m R N A translation. Slow-codons are represented by lattice sites w i t h a reduced particle movement rate p < p\. Placements of slow-codons or "defects" (blue shaded lattice sites), (a) Two adjacent defects subsumed in a L = 5 finite segment, (b) T w o defects separated by four normal codons and subsumed i n L = 11 finite segment. . 2 5.2 64 66 L i s t of F i g u r e s 5.3 5.4 x (a) Comparison of steady-state currents J\{p2) for a single defect derived from M C simulations w i t h those obtained from a L-segment F S M F method. For L > 5 , the F S M F results are w i t h i n 2 % of those from M C simulations. T h e boundary injection and extraction rates were not rate limiting (a = 8 = 1 0 ) . (b) Further reduction of the steady-state current as successive, identical defects are added. T h e first few defects cause most of the reduction i n current. Note that the ribosome current varies significantly as P2 is changed. Thus for clarity, the currents i n the figure have been normalized by the slow-codon hopping rate ( j ^ ) - Also note that for small P2, J2/P2 ~ 1 / 2 consistent w i t h E q . 5 . 3 6 8 (a) Steady-state currents across a lattice w i t h two identical interior defects spaced k sites apart. T h e current is suppressed most when the defects are closely spaced, (b) T h e density profiles near a pair of defects (p = 0 . 1 5 ) of various spacings k from F S M F computations. T h e circled points correspond to the defect positions for k = 6 . For larger k, density boundary layers decay, allowing particle build-up behind the second barrier, enhancing the current 7 2 (a) T h e dependence of the normalized steady-state current J2(P2, k)/Ji{.P2) as a function of the defect separation, (b). T h e upper bound for the mean current of a lattice w i t h m defects randomly distributed w i t h i n the central N = 3 0 0 sites. A lower bound as m — > N is J ( m — > N) AN —* P2/4 7 3 2 5.5 R 6.1 (a) Division of the a-hemolysin channel into the cis, vestibule, transmembrane, and trans regions. T h e D N A hairpin is inserted into the channel from the cis region. T h e number of D N A bases that can fit into the vestibule region is set to N = 1 6 , while the number of D N A bases than can fit into the transmembrane region is N = 1 2 . (b) A n example of the D N A hairpin translocation experiment. T h e process of D N A translocation is described by two coordinates n and ri\. n gives the position of the first hairpin base relative to the m o u t h of the transmembrane part of the channel, while index rii gives the number of hairpin base pairs that have been dehybridized. . . 7 9 T h e enumeration scheme for D N A base pairs i n the hairpin and the associated state space for the 2 D stochastic model. In this example N = 7 , NA — 4 , and there is a defect at base 4 . States w i t h n\ = 3 are inaccessible as the presence of a defect at site 4 makes it impossible to have exactly 3 base pairs dehybridized. Transitions between states w i t h different values of n have attempt frequency 00. In all cases, transitions between states w i t h different values of n have attempt frequency \x. States with n = 2 A f + A / l o o p + N (a D N A strand pulled completely into the trans chamber), and n = — N (a D N A strand pulled completely into the cis chamber) are absorbing 8 0 V L 0 0 6.2 x 0 0 L 0 A List of Figures 6.3 7.1 7.2 7.3 XI (a) T h e best fit predictions of the stochastic model when b o t h the defect free and defect-bearing strands are fitted simultaneously. A t voltages below the peak voltage, the D N A strands escape the pore p r i m a r i l y by sliding back into the cis chamber. A t higher voltages, the hairpin completely dehybridizes and slides through the channel into the trans chamber. T h e inset shows that the voltage corresponding to the mean first passage time peak can be controlled by altering the length of the hairpin's p o l y - A tail, (b) T h e probability of the D N A strands escaping into the cis chamber, (c) T h e probability of having UQ bases pulled into the transmembrane region at the time the hairpin fully dehybridizes (Piast)- Solid lines are the results for the defect-free 2-3 free energy surface, dashed lines are for the defectbearing 2-3 free energy surface (a) T h e geometry of the confining pore. Unless otherwise noted L = 1.5nm. (b) T h e wall potential as a function of radial position. T h e radius at which the water-wall potential energy is equal to fc^T is taken to equal the physical radius of the pore A d s o r p t i o n and desorption isotherms for water and non-ion-bearing pores. T h e solid lines correspond to the adsorption isotherms, while the dashed lines correspond to the desorption isotherms. Hysteresis was not observed in the 0.54 and 0.44nm pores. Lines are drawn as guides only. A t the density of liquid water at 298K (l.Og/cc) the expected number of waters adsorbed i n the 0.44nm pore is 8, 21 in the 0.54nm pore, 41 in the 0.64nm pore, and 67 i n the 0.74nm pore A d s o r p t i o n isotherms for water into ion-laden pores. T h e horizontal line gives the expected number of water molecules i n the physical volume of each pore at the ambient density of liquid water (m 1.0g/cm at 298K), while the vertical line gives the chemical potential of saturated water vapour at 298K ( « - 5 8 . 5 k J / m o l ) Characteristic water structures in pores w i t h geometric diameters of (a) 0.44nm, (b) 0.54nm, (c) 0.64nm and (d) 0.74nm. T h e line of water molecules in the 0.44nm pore is nearly perfectly hydrogen bonded w i t h an average of 1.8 hydrogen bonds per water molecule. \x\v = —56.0kJ/mol Water Oxygen - Water Oxygen radial distributions functions for pores containing pure water, and pores containing single ions (fiw = —56.0kJ/mol). T h e ions produce deviations from the pure water structure that are most pronounced i n pores w i t h smaller radii, and increase w i t h the size of the ion's charge. In the 0.74nm pore, the perturbation of the water structure resulting from the ion is small F i l l i n g i n the 0.44nm pore w i t h a C l ~ ion. (a) fiw = —66kJ/mol (n = 4) (b) fj, = — 6 4 k J / m o l (n = 6) (c) the low density state at nw = —62kJ/mol (n = 7) (d) the high density state at ji = —62kJ/mol (n = 11) (e) fi = — 5 9 k J / m o l (n = 12) (f) fi = —56kJ/mol (n = 13). A t fiw = —56kJ/mol states w i t h 12 water molecules similar to the state shown i n (e) and states w i t h 13 water molecules are equally probable. 90 102 Ill 3 7.4 7.5 7.6 113 115 117 w w w w w w w w w w 118 List of Figures 7.7 xii F i l l i n g i n the 0.44nm pore w i t h a K ion. (a) fiw = —66kJ/mol (nw = 5) (b) fjL = - 6 4 . 0 k J / m o l (nw = 7)(c) fiw = - 6 2 . 0 k J / m o l (nw = 9) (d) fi = - 6 1 . 0 k J / m o l (n = 10) (e) fi = - 5 9 . 0 k J / m o l (n = 13) (f) H = - 5 6 k J / m o l (n = 14) 119 7.8 F i l l i n g i n the 0.54nm pore w i t h a C P ion. (a) fiw = —66kJ/mol (nw = 6) (b) fiyv = —64kJ/mol (nw = 8) (c) the low density state at fiw = —61.1kJ/mol (nw = 15) (d) the high density state at fiw = —61.1kJ/mol (n = 31) (e) Lt = - 5 9 k J / m o l (n = 34) (f) fi = - 5 6 k J / m o l (n = 40). 120 7.9 F i l l i n g i n the 0.54nm pore w i t h a K ion. (a) Ltw = —66kJ/mol (nw = 7) (b) fiw = —64kJ/mol (nw = 9)(c) the low density state at fiw = —60.2kJ/mol (nw = 17) (d) the high density state at fiw = —60.2kJ/mol (nw = 26) (e) LL = —58.9kJ/mol (nw = 31) (f) JJL = —56kJ/mol (nw = 37) 121 7.10 F i l l i n g i n the 0.64nm pore w i t h a C l ~ ion. (a) Liw = —66kJ/mol (n = 6) (b) jj = —64kJ/mol (nw = 7) (c) the low density state at /J, — —60.6kJ/mol (nw = 19) (d) the high density state at fj,w = —60.6kJ/mol (n = 50) (e) fiw = —59.8kJ/mol (n = 53) (f) fiw = —56kJ/mol (nw = 63) 122 7.11 F i l l i n g i n the 0.64nm pore w i t h a K ion. (a) fiw = —66kJ/mol (nw = 6) (b) fiw = —61.1kJ/mol (nw — 14) (c) the low density state at fiw = —59.9kJ/mol (nw = 22) (d) the high density state at fiw — —59.9kJ/mol (nw = 49) (e) fi = —58.9kJ/mol (n = 52) (f) fi = —56kJ/mol + W W w W w w w w w w w w + W W w W w w w + w w w lnw = 60) 7.12 F i l l i n g i n the 0.74nm pore w i t h a C l ~ ion. (a) fiw = —66kJ/mol (nw — 5) (b) fi = - 6 4 k J / m o l (n = 6) (c) fi = - 6 0 . 3 k J / m o l (n = 20) (d) immediately prior to filling fi = —59.9kJ/mol (nw — 25) (e) immediately after filling fi — —59.8kJ/mol (n = 83) (f) fiw = —56kJ/mol (n = 93). Parts (a) and (b) show the tendency for water to layer i n the wall potential m i n i m u m even when only the ion's first solvation shell is present. w w w 123 w w w w w 124 7.13 F i l l i n g i n the 0.74nm tube w i t h a K ion. (a) fi — —66kJ/mol (n = 6) (b) fi = —64kJ/mol (nw = 7)(c) immediately prior to filling fi — —59.7kJ/mol (nw = 24) (d) immediately after filling fiw = —59.4kJ/mol (n = 82) (e) fiw = —57.5kJ/mol (nw = 88) (f) fiw = —56kJ/mol + w w w w w (nw = 92). In (f) a second inner layer of water is beginning to form 7.14 The total water-ion interaction energy (closed points), and the total particlewall interaction energy (open points) for the monovalent ions i n (a) 0.44nm, (b) 0.54nm, (c) 0.64nm, and (d) 0.74nm pores 7.15 The total water-water interaction energy (closed points), and the total configurational energy (open points) for the monovalent ions i n pores with radius (a) 0.44nm (b) 0.54nm (c) 0.64nm (d) 0.74nm 7.16 R a d i a l probability density profiles for C l ~ pores w i t h radii (a) 0.44nm (b) 0.54nm (c) 0.64nm (d) 0.74nm. In each case, note that the tendency for water to layer the inner surface of the pore is evident even at the lowest chemical potentials 125 126 126 127 List of Figures xiii 7.17 R a d i a l probability density profiles for K pores w i t h radii (a) 0.44nm (b) 0.54nm (c) 0.64nm (d) 0.74nm. Note that the filled state hydrogen distributions exhibit increased density at large radii compared to the profiles of + F i g . 7.16 7.18 R a d i a l distribution functions for K and C P i n pores w i t h radii equal to 0.44nm and 0.54nm. In b o t h cases \i\v — —59.0kJ/mol 128 + 129 7.19 (a) T h e radial pressures i n K (closed symbols) and C P (open symbols) bearing pores (b) the axial pressures. In b o t h cases lines are drawn as a guide to the eye. Errors are approximately 15% of the displayed values . . 133 7.20 T h e behaviour of the axial pressure PL under scaling of the pore length L. W h i l e the pressure i n the pure water pore oscillates as the length is scaled, the variation i n P over the range of displayed values of L is much smaller than the variation i n the K bearing pore ([xw = —57.0kJ/mol) 136 7.21 T h e grand potential of K (circles) and C P (squares) bearing pores with (a) R = 0.44nm, (b) R = 0.54nm, (c) R = 0.64nm, and (d) R = 0.74nm. Energies are i n ( k J / m o l ) of pore 138 7.22 T h e entropy i n K (circles) and C P (squares) bearing pores w i t h (a) R = 0.44nm, (b) R = 0.54nm, (c) R = 0.64nm, and (d) R = 0.74nm. Entropies are in ( k J / m o l • K ) of particles 139 7.23 (a) T h e adsorption isotherms for K bearing pores w i t h radius R = 0.54nm w i t h different lengths. T h e (b) ion, (c) water oxygen, and (d) water hydrogen radial probability profiles. T h e points are for the L — 1.5nm pore while the lines are for the L = 3.1nm pore (LIW — —58.2 k J / m o l ) 140 7.24 T h e (a) solute, (b) water oxygen, and (c) water hydrogen radial probability profiles for the different species in capped and periodic pores. K ° is a neutral solute w i t h the same Lennard-Jones parameters as K . In all cases the pore radius was R = 0.28nm and the pore was filled w i t h n = 34 water molecules 141 + L + + + + + w A.l A sample T A S E P and the heap structure used i n the B K L - b a s e d simulation. T h e value of each of the leaves equals the rate at which the associated particle move can take place. R e d shaded leaves w i t h a value of zero indicate moves that are impossible i n the current state of the lattice. Listed in the column on the right are the values of s, the random selector i n A l g . 6. T h e green hatched nodes are the nodes visited by A l g . 6 while selecting the next move, and the solid green node (Z ) is the selected move 4 159 xiv Preface T h e content of this thesis has been based, i n part, on the following articles: 1. "Totally asymmetric exclusion processes w i t h particles of arbitrary size", G . Lakatos, T . C h o u , J. Phys. A 36 (2003) 2027 (Chapter 3) T h e original idea for the research presented i n this paper was due to Professor T o m C h o u ( U C L A ) . T h e design of the research program was suggested by Professor Chou, while the methodologies and results were developed by G . Lakatos w i t h guidance from Professor C h o u . T h e final manuscript was written by G . Lakatos w i t h revisions and additions by Professor Chou. 2. "Steady-state properties of a totally asymmetric exclusion process w i t h periodic structure", G . Lakatos, T . Chou, and A . Kolomeisky Phys. Rev. E 71 (2005) 011103 (Chapter 4) T h e idea for the research presented in this paper was originated by Professor A n a t o l y Kolomeisky (Rice) and Professor T o m C h o u ( U C L A ) . T h e design of the research program, the methodologies, and the presented results were developed by G . Lakatos w i t h guidance and suggestions from Professors C h o u and Kolomeisky. T h e final manuscript was written by G . Lakatos w i t h revisions and additions from Professors C h o u and Kolomeisky. 3. "Clustered Bottlenecks i n m R N A Translation and P r o t e i n Synthesis" ,T. Chou, G . Lakatos Phys. Rev. L e t t . 93 (2004) 198101 (Chapter 5) The idea for the research presented i n this paper is due to Professor T o m Chou. Professor C h o u was also responsible for the design of the research program. The Preface xv methodologies were developed by G . Lakatos, the research and results are due to G . Lakatos w i t h guidance and suggestions from Professor C h o u . Professor Chou co-wrote the text w i t h G . Lakatos. 4. "First Passage times of driven D N A hairpin unzipping", G . Lakatos, T . Chou, B . Bergersen, and G . N . Patey Phys. B i o l . 2 (2005) 166 (Chapter 6) T h e research was proposed by G . Lakatos and Professors Bergersen and Patey. T h e design of the research program was developed by G . Lakatos and Professor Patey. T h e research and results are due to G . Lakatos. T h e manuscript was written by G . Lakatos w i t h suggestions from Professors Patey, C h o u , and Bergersen. 5. "Water A d s o r p t i o n in Ion Bearing Nanopores", G . Lakatos, G . N . Patey, To be submitted to the Journal of Chemical Physics (Chapter 7) T h e idea for the research presented i n this paper was conceived by G . Lakatos and Professor Patey. T h e design of the research program was developed by G . Lakatos and Professor Patey. T h e research and results are due to G . Lakatos w i t h guidance and suggestions from Professor Patey. T h e manuscript was written by G . Lakatos with suggestions from Professor Patey. As part of the b o d y of software used i n the preparation of this thesis, the following third-party numerical analysis codes were employed 1. T h e A r n o l d i m a t r i x diagonalization routines [1] employed i n Chapters 4, 5, and 6. T h e implementation used is included w i t h the commercial software package M A T L A B [2]. 2. T h e simplex multi-dimensional minimization scheme [3] employed i n Chapter 6. T h e implementation used is included w i t h the commercial software package M A T L A B . 3. T h e M u l t i p l e - F o l d ( M F O L D ) 2-state hybridization server used i n Chapter 6 [4]. Preface xvi A l l other simulation and analysis codes, including the M o n t e Carlo simulation programs used in Chapters 3, 4, 5 and 7, were developed by G . Lakatos. T h i s includes extensive code and scripting developed to drive the third-party numerical analysis codes enumerated above. XVII Acknowledgements W r i t i n g a thesis is a lengthy process that cannot be completed without support and guidance from many individuals. I would like to begin by thanking my supervisors, Professors Patey and Bergersen. T h e y have patiently guided and supported my research throughout my time at U B C , and provided valued encouragement when difficulties arose. I would also like to thank my collaborators, and i n particular Professor T o m Chou, who has provided helpful suggestions and good advice. T h e University itself has provided a wonderful environment i n which to learn about the process of doing research, and through the University Graduate Fellowship program, has provided extensive financial support. Finally, I would like to thank the members of my family for their care and encouragement. 1 Chapter 1 Introduction Over the last decade there has been an explosion of interest i n molecular transport processes and structure formation on extremely small length scales. M u c h of this interest has been driven by the rapidly advancing fields of nano- and biotechnology. In this thesis, several models of structure formation and particle transport applicable to nanoscale channels are developed, and investigated using analytical methods and computer simulations. We first focus on transport in the biological domain where there are several transport processes that are fundamental to the functioning of a l l living organisms. Examples include the movement of so-called motor proteins such as myosin [5] that is ultimately responsible for all muscular movement i n multi-celled life. M y o s i n molecules, along w i t h many other motor proteins, have been found to move i n one direction along long polymerlike fibres [6]. B y changing their molecular conformation, myosin molecules are able to make discrete and equally-sized steps from one binding site to the next along the length of the fibre. In addition to the various motor proteins, we can consider the processes of gene expression, and protein synthesis [7]. Here we find multiple distinct transport processes that are critical elements i n the mechanism that converts the information stored in an organism's D N A into functional proteins. A s an example, consider the final step of protein synthesis where a molecule of so-called messenger R N A is used as a template to construct a protein [8]. T h i s final stage of the protein production process is performed by ribosome molecules which bind to one specific end of the polymer-like messenger R N A and then, through a series of discrete and equally sized steps, proceed to slide along the R N A until they reach the opposite end of the molecule. W i t h each step a ribosome makes, an amino acid is added to the growing protein that is attached to the ribosome. T h e movements of motor proteins and ribosomes have been investigated using a variety of theoretical and experimental techniques. For example, first principles quantum Chapter 1. Introduction 2 calculations have been used to study the interactions between myosin and the adenosinetriphosphate molecule [9], as well as the binding of amino-acid carrier molecules ( t R N A s ) to ribosomes. Similarly, molecular dynamics simulations [10] have been used to study the conformation changes involved i n the stepping behaviour of motor proteins. While these single molecule studies have proven very useful in understanding the behaviour of individual motor proteins or ribosomes, they have had less success modeling the collective behaviour of the multiple molecules that are located on the transport pathway at any one time. Moreover, these studies tend to be highly specific to particular molecular transport systems, and are difficult to generalize to other biological transport processes. T h i s latter point is particularly acute given that many biological transport processes share a number of pertinent features. For example, the movement of b o t h ribosomes and myosin molecules occurs along stand-like filaments and so can be treated as one-dimensional. Additionally, i n b o t h cases the molecules only move i n one direction along the filaments. Finally, myosin molecules and ribosomes both move i n a series of discrete steps, all equal in length. A model which captures these three salient features is the so-called Totally Asymmetric Simple Exclusion Process ( T A S E P ) [11-35]. Originally developed to model the movement of ribosomes [11], the T A S E P has been successfully applied i n a variety of contexts ranging from the study of highway traffic [19, 20, 29], to the foraging behaviour of ants [35]. A s the T A S E P model forms the basis of the first part of this thesis, Chapter 2 contains a description of the model and outlines the techniques used to analyze it. Chapters 3 and 4 describe a set of modifications to the T A S E P designed to enhance its utility in modeling biological transport processes. T h e extensions of the T A S E P and the solution techniques developed in Chapters 3 and 4, are then applied i n Chapter 5 to study the movement of ribosomes along messenger R N A ( m R N A ) molecules where the movement rate of the ribosome varies w i t h its position along the length of the m R N A . Another example of transport i n a biologically derived system is the translocation of D N A and R N A through membrane embedded nanopores. T h e translocation of D N A and R N A through membrane pores occurs naturally i n all organisms containing a cell nucleus, and is a critical step in the viral parasitism of cells [36]. Recently, the process of D N A translocation has been the subject of extensive experimental investigation [37-41], Chapter 1. Introduction 3 which has established that the characteristic time scale for the sliding of D N A strands with approximately one hundred bases is on the order of milliseconds. Theoretical work [42-47] has resulted i n the development of several models of the thermodynamics of the translocation process. T h i s includes descriptions of how the structural properties of single stranded D N A , and the interactions between the D N A and the transmembrane pore, affect the rate of translocation. A separate line of investigation [48-51] has considered the process of dehybridizing double-stranded D N A and related the force required to separate two stands of hybridized D N A to the sequence of the D N A . C o m b i n i n g the work done on D N A translocation and dehybridization, recent experiments [52, 53] have investigated the process of inducing D N A dehybridization by forcing the D N A through narrow channels. This process of driven D N A dehybridization has the potential to become the physical basis of novel D N A sequencing devices. In Chapter 6, a stochastic model is developed that couples the D N A dehybridization and translocation processes. T h e model is used to investigate the variations i n the translocation time of a sample set of double stranded D N A molecules as the voltage driving them through a transmembrane pore is varied. Structure formation in nanoscale channels is investigated i n Chapter 7, where a study of water adsorption i n narrow cylindrical pores is presented. T h e behaviour of water in the interior of narrow, ion-bearing pores is of tremendous interest from the biological perspective due to the importance of transmembrane ion channels i n maintaining cellular function [54-59]. Molecular dynamics studies have, for example, found that the ordering of water in the narrow confines of cylindrical pores produces a decrease i n the ion mobility by a factor of two to four compared to the bulk mobilities [60]. Additionally, there has been significant recent interest i n utilizing carbon nanotubes as fluid storage devices and conduits [61-64], and here the structure of confined water can influence the carrying capacity of these nano-containers. In Chapter 7, we use M o n t e Carlo simulations in the grand canonical ensemble to investigate the effects of adding a single ion to the interior of a nanometer scale cylindrical pore on the adsorption of water. Finally, i n the last chapter of the thesis we review and summarize our results. A d ditionally, we comment on many of the outstanding questions posed by the work in the thesis, and future research directions based on the results of the thesis are proposed. 4 Chapter 2 The Totally Asymmetric Simple Exclusion Process The totally asymmetric simple exclusion process ( T A S E P ) is a lattice gas model that describes the movement of particles along a one-dimensional pathway. The standard T A S E P consists of a one-dimensional lattice of N sites w i t h open boundaries at either end (see F i g . 2.1(a)). Particles exactly one lattice spacing i n size can enter the lattice at the left end w i t h an attempt rate a, provided that the first site of the lattice is unoccupied. Once in the lattice, particles move to the right in steps exactly one lattice site in length, and make these steps w i t h an attempt frequency p . A s the particles have a hard core interaction, no two particles can simultaneously occupy the same site. Thus, a particle can only move right if the site immediately to the right of its present location is unoccupied. A n y particle in the last site of the lattice is removed w i t h an attempt rate B. In the variant of the T A S E P discussed i n this thesis, multiple particle movement events are forbidden from occurring at exactly the same time. However, the time separating two successive events, while finite, may be arbitrarily small. T h i s restriction on the movement, injection, and extraction events, produces the so-called r a n d o m - s e q u e n t i a l update procedure [65]. A s used in this thesis, the T A S E P has three noteworthy features; it is a strictly one dimensional model, it is discrete in b o t h the transported objects (the particles) and the steps they take, and it's unidirectional as all the particles move left to right. A s will be shown in this first part of the thesis, these properties make the T A S E P a very valuable tool for modeling biophysical transport processes. In fact, early investigations [11, 12] of the T A S E P focused on using the T A S E P as a model for the movement of ribosome molecules along strands of messenger R N A . In this context, the primary quantities of interest were the steady-state particle current J, defined as the number of particles passing an arbitrary point i n the lattice per unit time, and ctj, the particle density at any particular lattice site i. T h e density at a site i was defined to be the probability of finding the site occupied Chapter by a particle. 2. T h e T o t a l l y A s y m m e t r i c S i m p l e E x c l u s i o n Process 5 These early studies used approximate analytic methods, and computer simulations, to determine the steady-state behaviour the T A S E P under different values of the entry and exit rates, given a fixed value of the internal hopping rate. W h i l e the early investigations of the model d i d not produce an exact solution, the T A S E P was shown to exist i n one of three solution regimes, or phases, depending on the relative values of the entry and exit rates (see F i g . 2.1(b)). In the case where a < ft, and a < p/2 the T A S E P was i n a low density entry-limited regime. In this case, getting particles into the lattice was the step that limited the overall rate of particle transport. Similarly, when P < a and P < p/2, the T A S E P was i n a high density exit-limited regime where removing particles F i n a l l y when a and P were both at the right end of the lattice limited the current. greater than p/2 the T A S E P was i n the maximal current regime. In this case, the internal motions of the particles rather than injection and extraction at the boundaries, limited the particle current through the lattice. T h e boundaries between the phases shown in F i g . 2.1(b) are sharply marked by discontinuities i n measurable properties of the T A S E P . Specifically, at the transition between the entry and exit limited phases (the a = P line), the T A S E P current ( J ) is continuous, while dJ/da, dJ/dp, and the mean bulk density of the T A S E P lattice a l l show a j u m p discontinuity. Similarly, at the boundary limited to m a x i m a l current transitions (the a = p/2 and P = p/2 lines), the bulk density and current are continuous, while the first derivative of the density and the second derivative of the current show discontinuities. B y convention, the transitions between the solution regimes are classified by the order of the discontinuity i n the current and thus the entry to exit limited transition is said to be first order, while the boundary limited to maximal current transitions are second order. After the early investigations of the T A S E P , the next major advance i n the study of the T A S E P was the development of exact solutions by D e r r i d a et al. [15] and Schiitz and Domany [16]. T h e m a t r i x product solution of D e r r i d a i n particular provides concise expressions for the mean densities, density correlations, and currents i n the steady-state of the T A S E P . For example, the current i n a Af-site T A S E P is given by the expression j RN-I(p/P) = ~ P - RN-iip/a) R (p/P)-R (p/a) N N (2-1) Chapter 2. The Totally Asymmetric Simple Exclusion T T T T T Process 0 i (b) 1 ! J V 1 1 Exit Limited Phase 0.6 Figure 2.1: (a) T h e standard totally asymmetric simple exclusion process ( T A S E P ) . Particles are injected on the left at a rate a, traverse the lattice from left to right by making hops w i t h an attempt rate p , and are removed from the lattice w i t h an attempt rate B. N o more than one particle can occupy a lattice site at a time, (b) T h e three solution regimes of the standard T A S E P . Note that here p has been set equal to 1. I N S E T S : T h e density profile along the T A S E P lattice. In the entry and exit limited phases two density profiles are possible depending on the values of (a,/?) relative to the a = 1 — B line (the dashed line i n the plot). Chapter 2. The Totally Asymmetric Simple Exclusion Process 7 where, R M ) - " T H (x)-2^ N U - 1 ) , m { 2 N N + -_ i 1 ) j l ) l x* <- ((2- )2) 2 2 In the limit of large JV, the exact expression for the current takes on one of the three phases shown i n F i g . 2.1(b) and shows the same discontinuities at the phase boundaries as were known to M a c D o n a l d et al. [11, 12]. Additionally, when N is large the matrix product solution generates the following simple expressions for the current and the density in the interior of the T A S E P lattice i n the three solution regimes (Entry) a < p/2 J = a(l - a/p), a = a/p J = p(l-p/p),o = L L a<P (Exit) 0<p/2 R R (2-3) 8 < a (Maximal) l-p/p a > a* = p/2 J m a x = p/4, a = 1/2 m a x The expressions for the current i n E q . 2.3 were known long before the exact solution to the T A S E P was found [11, 12]. In fact, simple assumptions lead very readily to these results. Consider the entry-limited case; we assume the density is approximately constant throughout the T A S E P lattice and label this density OL- W e then set OL equal to the ratio of the injection rate to the internal stepping rate, ai — ct/p. T h e particle current into the first site i n the lattice is then J L JR. = P(l — P/p) i n the exit-limited phase. = a(l — a/p). A similar argument gives In the m a x i m a l current phase we assume that the density is constant and the probability of any particular site i n the lattice being occupied is uncorrelated w i t h the probability of an adjacent site being occupied. If we let the density be a , then the assumption of no correlation between adjacent lattice sites gives J m a x = pa (I — a). M a x i m i z i n g this expression gives a = 1/2 and J = p/A. A s we ignore correlations between lattice sites when producing these results, these simple calculations are so-called mean field calculations. Chapter 2. T h e Totally Asymmetric Simple Exclusion Process 8 Another important feature of the T A S E P that manifests itself i n E q . 2.3 is the symmetry between the entry- and exit-limited phases. Consider the entry-limited solution in Eq. 2.3 and replace a w i t h (3. T h i s immediately returns the correct current for the exit limited phase. Similarly, the exit-limited density equals 1 — o _ , under the substitution a —> (3. T h i s is a consequence of so-called particle-hole symmetry in the T A S E P , where the movement of a particle from left to right in the entry-limited phase is the same as the movement of an empty site, or hole, from right to left i n the exit-limited phase. To see why a seemingly simple model like the T A S E P is so difficult to solve, consider the equations for the mean occupancy of a lattice site. Let Xi(t) be the state of lattice site i at time t, where Xi(t) = 1 if the lattice site is occupied, and Xi(t) = 0 if the lattice is unoccupied. T h e mean occupation of an interior lattice site and the first lattice site are then given by d(xi) = p((Xj-l) - (Xi) + {XiX ) - i+l (Xi-iXi)) (2.4) = a(l-Ei)-p(xi(l-x )>. 2 Here, we can note that the symmetry-exchange operation interchanges a and f3 and replaces Xi w i t h (1 — x -. i) N dt d t = = i n E q . 2.4, yielding i+ P((x -i) ~ N ( Z j v - i + l ) - {x -iX - i) N + N i+ (x _ iX _ )) N i+ N i+2 , (2.5) -p(x )+p{x -i(l-x )). N N These are the correct equations for an interior lattice site N xjv-i+i and the last site in the lattice. Note however, that while Eqs. 2.4 and 2.5 are equations for the mean occupancy of single lattice sites, they depend on the occupancy correlations for adjacent lattice sites. Similar equations can be written for the two-site occupancy correlations, however these equations include three-site correlation terms. T h i s pattern continues for the higher correlations and so ultimately knowledge of the Af-site correlation functions is needed in order to produce an exact solution, or correlations of some order must be ignored resulting Chapter 2. The Totally Asymmetric Simple Exclusion Process 9 in a mean field solution. Unfortunately, while the mean field results provide exact expressions for the currents and the densities of the T A S E P lattice far from the boundaries, mean field approaches fail i n capturing the exact density profile along the T A S E P lattice (insets in F i g 2.1(b)). This can be explained if one assumes that the inter-site correlations i n the T A S E P lattice decay to zero far from the lattice boundaries. If this were true, then i n the interior of the lattice mean field theories would be exact. A s there are no particle sources or sinks in the interior of the T A S E P lattice, the particle current must be the same everywhere, and thus calculating the current i n the region where mean field theory is exact will yield the one correct current. Conversely, in the regions near the lattice boundaries where there are non-zero inter-site correlations, mean field methods would not exactly capture the variation in the particle density as the boundaries are approached. In fact, this is the case and the correlations between adjacent lattice sites i n the T A S E P are found to decay to zero as the inverse of the distance from the boundaries i n the m a x i m a l current phase, and exponentially i n the boundary limited phases [15]. T h e deviations from the bulk density near the boundaries (the boundary layers) also decay to zero, decreasing as the inverse of the square root of the distance from the boundaries in the m a x i m a l current phase, and exponentially in the boundary limited phases [15]. A s mentioned above, while mean field approaches have been useful i n studying the standard T A S E P , they do not successfully capture all the interesting aspects of the model. T h i s problem is particularly acute when dealing w i t h modified variants of the T A S E P for which no exact solution exists, and in these cases having a benchmark for the quality of the mean field estimates is advantageous. In such situations, the standard means of studying the T A S E P and T A S E P - l i k e models has been to employ computer simulations based on the Monte Carlo method [66]. Given a system that can be found in one of a set of states S, the M o n t e Carlo method is a means of accurately calculating the probability of finding the system in any particular state St. F r o m a knowledge of the state probabilities, any property of the system can i n principle be computed. To apply this method to the T A S E P , first note that the state of the T A S E P is completely specified if we the know the occupancy of each lattice site X\ through x^. E a c h possible configuration of the T A S E P Chapter 2. The Totally Asymmetric Simple Exclusion Process 10 lattice ( x i , X 2 , x s , . . . , £jv) can be labelled w i t h a single index k, and then we have the following equation for the time evolution of the probability distribution for the lattice state ,N (2.6) In E q . 2.6, is the rate at which a T A S E P in state Si changes to state Sk- For example, if states Si and Sk differ by the movement of a particle i n the interior of the lattice, then Tik = p . E q u a t i o n 2.6 is the so-called master equation for the T A S E P . In this thesis we do not study the time evolution of the T A S E P , and focus exclusively on steady-state properties of the model. To compute the steady state probability of being in state k, the left hand side of E q . 2.6 is set to zero. M a n y different M o n t e Carlo algorithms for determining the steady-state probabilities of the lattice states Si are available, and i n this thesis we have used a continuous-time Monte Carlo algorithm developed by Bortz, Kalos and Lebowitz ( B K L ) [67]. Given a T A S E P lattice i n state S i , the B K L algorithm proceeds as follows 1. Compute the total transition rate (R) out of state Sk, where R = 2~2i=oi^k kir state Si associated w i t h a non-zero Any is accessible from state Sk i n a single particle movement. 2. Select and perform the transition to state Si w i t h probability T k i / R 3. Increment the total time spent in state Sk by 1/R and return to Step 1. A n efficient implementation of this algorithm for the T A S E P is described i n A p p e n d i x A. After a large number of iterations of the algorithm, the steady state probability of finding the T A S E P lattice in any state Sk will be well approximated by the ratio of the time spent in state Sk to the total simulation time. Note that i n the B K L algorithm, the probability of making a particular transition is only a function of the current state of the T A S E P lattice; previous transitions have no effect on this rate. Also note that as a consequence of the way the T A S E P was defined, given any two lattice states Si and Sk there is a finite sequence of transitions that will evolve the lattice from one state to the Chapter other. 2. The Totally Asymmetric Simple Exclusion Process 11 Given these properties of the T A S E P model and the B K L algorithm, it can be shown [68] that after a large number of iterations the probability distribution generated by the B K L algorithm is guaranteed to satisfy the steady state form of E q . 2.6. After each iteration of the B K L algorithm, statistics for the current J and density o~i can be taken. For the current, the number of particles that reach the last site of the lattice and subsequently leave the lattice is recorded. D i v i d i n g the total number of particles that exit the lattice by the total simulation time yields the M o n t e C a r l o approximation to the steady state current. Similarly, the densities can be computed by recording the state of the lattice at regular intervals of simulation time. T h e density at lattice site i is then approximated by dividing the number of samples'where the lattice site was found occupied by the total number of samples taken. In this first part of the thesis, mean field methods and M o n t e Carlo simulations are used to investigate the steady state currents and densities i n two modified T A S E P models. Chapter 3 removes the size l i m i t a t i o n the standard T A S E P applies to the translating particles, while Chapter 4 extends the T A S E P to the case where the hopping rates vary along the length of the T A S E P lattice. In the spirit of M a c D o n a l d ' s original work on the T A S E P , the methodologies developed to treat the variable hopping extension of the T A S E P are utilized i n Chapter 5 to explore the effects of variable codon usage on the translation rate of messenger R N A . 12 Chapter 3 Totally Asymmetric Exclusion Processes with Finite Size Particles 3.1 1 Background As described i n Chapter 2, the standard T A S E P consists of particles that make steps with a length equal to the particle diameter. However, i n numerous applications, the size of the particles and the size of the steps the particles make are not identical. In other words, the particle may be large and occlude some number d > 1 lattice sites but only take steps of one lattice site. Specific biophysical examples include ribosomes moving along m R N A and large molecules or vesicles that are shuttled by motor proteins along microtubules. Ribosomes are approximately 20nm i n diameter, but move one codon at a time taking steps of approximately 2nm. Kinesin motor proteins that carry vesicles are typically greater than 50nm in size, and move along microtubules w i t h steps that are roughly 5nm long. In each of these examples there is a driven system of particles that take discrete steps along a one-dimensional track, and hence it is natural to attempt to model the particle motions using a T A S E P - l i k e model where d > 1. Steady-state properties of the d — 1 T A S E P such as currents, densities, and their correlations have been computed exactly using matrix product and recursive methods. However, when the size of the particles is greater than the size of a single lattice site, the matrix product method has yet to be successfully applied. In this chapter, we use mean field methods to derive approximations for steady state currents and densities for a T A S E P w i t h d > 1. O u r approach is guided by the results of M o n t e Carlo simulations that indicate the qualitative features of the T A S E P phase diagram remain intact when d > 1. Thus we assume that the steady-state currents exist i n one of three regimes, ^ h e work in this chapter is based on "Totally asymmetric exclusion processes with particles of arbitrary size", Greg Lakatos, Tom Chou, J. Phys. A 36 (2003) 2027-2041 Chapter 3. Totally Asymmetric Exclusion Processes with Finite Size Particles 13 exit-limited (with currents determined by /3), entry-limited (with currents determined by a ), and the m a x i m a l current phase (with currents determined by the interior hopping rate). We first consider simple mean field theories that do not explicitly take into account the hard-core exclusion of the large particles. We then consider refined mean field methods where the hard-core exclusion is taken into consideration by explicitly enumerating the possible states of the T A S E P lattice. These mean field results are then compared to results from Monte Carlo simulations. 3.2 M e a n Field Theories Consider a T A S E P of identical particles with diameter d > 1 lattice sites that are driven through a one-dimensional lattice w i t h N ^> d sites. T h e particles interact through hard core repulsion. T h e entry rate a, and exit rate (3 are all normalized by the internal hopping rate p. A s a result, all current results will be in units of the internal hopping rate. In all the work that follows, the location of a particle's left edge is used to determine the location of the particle in the T A S E P lattice. 3.2.1 Simple Mean Field Methods The simplest mean field approximation to the steady-states of this system assumes that the particles are on average uniformly distributed along the lattice without correlations arising from their hard core repulsive interactions. T h e left edge of a particle of size d can move forward only if the d sites ahead of it are unoccupied. B y ignoring the correlations we can immediately write: J = Oi(l - Oi+i)(l - c r i + 2 ) • • • (1 - o i + d ). (3.1) Here <7; is the probability finding the left edge of a particle at lattice site i. Note that this simple mean field theory ignores the perfectly correlated exclusion between particles Chapter 3. T o t a l l y A s y m m e t r i c Exclusion Processes w i t h F i n i t e Size Particles 14 w i t h sizes greater than a single lattice site; namely that when a particle is at site i, the next (d — 1) sites are empty w i t h probability one. We now consider the case where the particle is in the interior of the lattice far from the boundaries. Guided by the results for the d = 1 T A S E P , we assume the density in this region is constant and set <7; = o for all i, giving J = o(l-o) . (3.2) d The m a x i m u m current, and its associated density is then found by setting d J / d a = 0. This yields d 1 d J m a x = ( d + l)^+i' °" m a x = (d + T) • ( 3 , 3 ) Now consider particle entry at the left end of the lattice, or particle exit from the right Figure 3.1 shows two possibilities for the rules of entry and exit when d > 1. For end. example, a particle may only enter incrementally, one lattice site at a time (Fig. 3.1(a)), or it may enter completely i n one step when the first d sites of the lattice are empty (Fig. 3.1(b)). However, i n the large limit, the steady-state current and densities are independent of the different entrance/exit schemes. To see this, note that partial entry is equivalent to complete entry into a lattice that is extended by d — 1 sites at the entrance boundary, while partial exit is equivalent to complete exit i n a lattice extended by d — 1 sites at the right boundary. Since d / N —> 0, the differences arising from the variations in entrance and exit are finite-size effects that vanish i n the limit of large N. For the sake of simplicity, we w i l l assume complete entry (Fig. 3.1(b)) i n all of our subsequent work. In the simple mean field method we will consider complete exit ( F i g . 3.1(d)), while in the refined mean field method described later partial exit ( F i g . 3.1(c)) will be used. To find the entry-limited current i n the simple mean field approximation, we enforce current continuity at the left edge of the lattice giving d d aY[(l-a ) j 3=1 = a ]l(l-a ). 1 1+j (3.4) j =l Setting Oi = o~\ = a for all i, the occupation probability near the left entry site is then Chapter 3. T o t a l l y A s y m m e t r i c Exclusion Processes w i t h F i n i t e Size (a) Particles 15 (c) ! I I ! i 1 i 1 1 ; ! ! I ! .... .• -Y i i i '• i i (b) (d) i rr-l r—r—r-1 r - m i-rnri m n m a ~ P Figure 3.1: E n t r y and exit mechanisms, (a) Partial entry (b) complete entry into the first d sites (c) partial exit (d) complete exit as the particle reaches the last site. o = o L = a and the entry-limited steady-state current becomes J The = a(l-a) . (3.5) d L T A S E P will be i n the entry limited phase with current Jr, provided f3 is large enough to ensure that extraction is not the rate-limiting step in the movement of particles from one end of the lattice to the other. We also require that J L < J m a x , giving an estimate for a* = l/(d + 1). Here a* is the boundary between the entry rate-limited and maximal current regimes. W h e n the current is exit rate limited we again enforce current continuity yielding: d-l P&N-d+l J_(l - d &N-d+l+j) = ON-d i=l Jj£(l - CTN-d+j)- (3.6) j=i As before we set ctj = a for all i. This gives J = P (l-(3). Setting E q . 3.7 equal to E q . 3.3 gives 0* = d/(d+ The can (3.7) d R 1). steady-state currents and densities produced by the simple mean field approach be summarized as Chapter 3. (Entry) T o t a l l y A s y m m e t r i c E x c l u s i o n Processes a < -^-j- a(l-a) d (Exit) ( 5 < J w i t h F i n i t e Size = a(l - a) o = P (l-P) a = l-P d L Particles = a L <B (l- B) d ^l d J R R P {l - P) < a(l - a) d (Maximal) 16 (-) d a > a* = — 1 3 d d j _= ^ = _ 1 8 _ d + 1" The currents, boundary densities, and interior densities i n E q . 3.8 reproduce the known results for the standard T A S E P when d = 1. However, as the densities are assumed constant this mean field method does not produce accurate density profiles. Additionally, the boundary densities are expected to be inaccurate since they are independent of d. The predicted currents (Eq. 3.8) are continuous across the (a, P) parameter space, while exhibiting a discontinuity in the first derivative at the transition between the entry and exit limited regimes, and discontinuities in the second derivative at the a = a* and P = P* lines dr dr 2 da , 2 da 2 2 = A y=a(*)- da 2 v = r v ( * ) dP\ d~ d /3=/3(*) id- l) ' ' d 2 (3.9) In later sections of this chapter the currents, densities, and phase boundaries produced by this simple mean field approach will be compared w i t h more accurate mean field theories and numerical simulations. Chapter 3. T o t a l l y A s y m m e t r i c E x c l u s i o n Processes w i t h F i n i t e Size Particles 17 Figure 3.2: (a) A typical configuration of particles associated w i t h the calculation of Z(n, L) given by E q . 3.11. T h e left edge of the particles is noted by a black dot, while the exclusion zone associated w i t h each particle is shown w i t h a blue rectangle (b) A typical configuration for the calculation of Z(n. L — 1). 3.2.2 Refined Mean Field Theory for the Maximal Current Phase We will see that the simple mean field approximation developed in the previous section is inaccurate for d > 1. A more refined approximation for a T A S E P w i t h locally uniform particle distributions is to consider a means of explicitly including the extended exclusion due to the length of the particles. A g a i n marking a particle's location by its left edge, consider a particle located in the lattice interior at site j, and consider an L-site section of the T A S E P lattice located immediately to the right of the right edge of the particle (see F i g . 3.2). Here we take L » d and L <C N. T h e probability that the particle is free to move to the right is equal to the probability that site j + d is empty. This can be calculated from the ratio P{T l+l > J + d\Ti = j) = Z{n,LZ(n, L) 1) (3.10) ' In E q . 3.10, n is the position of the left edge of particle i and Z(n,L) of ways of arranging n particles of size d in a lattice of L sites. is the number Note that Z(n, L) is Chapter 3. T o t a l l y A s y m m e t r i c E x c l u s i o n Processes w i t h F i n i t e Size Particles 18 exactly the canonical partition function for a set of n size d hard particles on the L-site sub-lattice. Z(n, L) can be written Z ( n , L ) = ( £ - " n - 1 ) ) n (3.11) which gives P(r l + 1 > j + d\n = j) Z(n,L- = 1) Z(n,L) L — (d — l ) n — n 1— L-(d-l)n ad (3.12) l-a(d-l)' In E q . 3.12 we have assumed the number density of particles i n the lattice interior is approximately constant and set a = n/L. J = P( n The maximal current J m a x = j)P( > n+1 T h e steady-state current is then 3 + d\ = j) = <^ ^~^ _iy (- ) d n 3 1 and its associated bulk density a m a x 13 is found by maximizing J with respect to a yielding 1 a Vd(Vd+l) : n (3-14) x Jry (Vd+1) ' 2 A s expected, E q . 3.14 reduces to the standard T A S E P result [13, 15] when d — I. Also i n the limit of large d, J finding the m a x and a left edge of a particle at m a x ~ 1/d. Note that a m a x is the probability of any particular lattice site i n the interior of the lattice. However, any particular lattice site i will be occluded by some part of a particle whenever the left edge of a particle is located at a site j where (i — d) < j < i. Thus under the assumption of constant cr, the probability that an interior lattice site is occluded by some part of a particle is da . max T h e steady-state current will be the m a x i m a l current, J m a x , as long as the injection and extraction rates of an open boundary T A S E P are large enough to ensure that internal moves are the overall rate limiting steps. A s w i l l be demonstrated later in the chapter, the current and density derived from this refined mean field approach Chapter 3. T o t a l l y A s y m m e t r i c E x c l u s i o n Processes Particles 19 i =i (a) a w i t h F i n i t e Size J l 0 0 0 0 h 1 0 0 0 3 0 1 0 0 0 4 0 0 1 0 0 0 5 0 0 0 1 0 0 s S S S (b) i = N-4 p 0 N-3 N-2 N - l N 0 0 1 s 0 0 0 S ! 1 0 0 0 S i 0 1 0 0 S 1 0 1 0 S i 0 J_ 0 l 0 IP 2 3 4 5 Figure 3.3: State enumerations for particles of size d = 4 near the lattice boundaries, (a) T h e d+1 distinct states at the entrance region, (b) T h e d+l distinct states of the exit region. Note that a 0 indicates the absence of the left edge of a particle at the lattice site. Thus in the states displayed i n (b), all the 0 sites may be occluded by part of a particle, but cannot carry a left edge. In both (a) and (b), all the states that are w i t h i n d moves of state s\ are included. match extremely well w i t h those found from Monte Carlo simulations run in the maximal current phase. 3.2.3 Refined Mean Field Theory for Boundary-limited Currents Here, we develop a mean field theory that accurately models the behaviour of the particles at the ends of the lattice. First consider the case where a is small and examine the entrance site at the left of the lattice as shown in F i g . 3.3(a). A s we assume that a size d particle immediately occupies the first d sites of the lattice after injection, enforcing current continuity in the first d lattice sites gives Chapter 3. T o t a l l y A s y m m e t r i c E x c l u s i o n Processes = aP{n >d) J L = Pin = w i t h F i n i t e Size Particles d,r > 2d). 20 (3.15) 2 Recall that the simple mean field assumption used to derive Ji and JR (Eqs. 3.5 and 3.7) completely ignores exclusion beyond one lattice site. Here we correct this deficiency by defining a relative probability (9) that a site w i t h i n the first d sites of the left edge of the T A S E P lattice is occupied. Note that the extended size of the particles then limits the number of possible states of the first d sites of the lattice to exactly d + 1 (states Si through s 5 for d — 4 particles shown i n F i g . 3.3). To approximately account for the effects of the extended exclusion of the d > 1 particles, we normalize the relative occupancy probabilities by the sum of the relative probabilities of each of the d+ 1 occupancy states. This gives The normalization factor d9(l — 9) ~ d + (1 — 9) x d is the sum of the probabilities of all possible states of the first d lattice sites (Fig. 3.3(a)). T h e probability that the first particle is at site d, and is not prevented from moving one lattice site to the right by the presence of another particle is P We ( r i = d , T 2 > 2 d ) = ^ ( i _ ^ - ^ i _ f l ) j . have assumed that 9 is approximately constant over the first d sites. (3.17) The mean probability of occupation of the first d sites is then 9(l-9) ~ d ° L ~~ de(i - ey- 1 9 l + (i - ey ~ i + (d-i)e' ( 3 ' 1 8 ) Using E q . 3.15 we find that 9 = a, thus in the entry limited phase a l)a a(l — a) 1 + (d-l)a 1 + (d- Jr, = (3.19) In E q . 3.19, oi is the mean probability of finding site i occupied by the left edge of a Chapter 3. T o t a l l y A s y m m e t r i c E x c l u s i o n Processes w i t h F i n i t e Size Particles 21 particle where 1 < i < d. A n approximate expression for the probability of finding site i < d occluded by some part of a particle is then For io . L the exit-limited case, we consider all possible states of sites N — d+1 through N. Note that for small values of (3, we expect particles to be close packed at the right edge of the lattice. A s a result, the left-edge occupation densities w i l l split into two classes. Lattice sites w i t h index i = N — kd for integer k > 0, are the sites where the left edge of the particles spend the majority of the time i n the close packed state. A s a result we expect the occupation densities at these sites to be higher than at the intervening sites by roughly a factor of 1/(3. To proceed, we assume the density variation of each of these two classes of sites is small near the exit boundary. Referring to F i g . 3.3(b), we then apply current continuity to transitions from state s to state S\. T h i s gives the following 5 relation for the steady-state current J R = (3P(r n = N) = P{r n = N - 1). (3.20) Since in each block of d sites there is only one site of the high occupancy class, we set P(xyv-fcd) — ON — OR without any additional normalization factor. Now, the probability that the last particle is at l — 1 is " ' (i - ey-i + (d - i)9(i -ey-*' ^ - n > where 9 is a mean relative probability of finding the left edge of the particle in one of the low-occupancy class sites. T h e quantity ( l _ 0 ) d - i + ( _ _ i ) 0 ( i -e) ~ , d (3.22) 2 is then the partition function for the last set of low-occupancy sites N — d+l,N 2,...,N-1. N relationship between 9 and 9 , which gives 6^(1 - 9 )(1 N - O)^ N To find another we can equate the current into and out of the state s N find = 9/(1 + (d-2)9). Solving E q . (3.20), we find (36 /{1 -6 ) — d+ 1 = 9(1 - 9) ~ (l d 2 - 9 ). N 3 Solving these two equations we Chapter 3. T o t a l l y A s y m m e t r i c E x c l u s i o n Processes w i t h F i n i t e Size Particles 22 1-/3 OR l + {d -l)p' JR (3.23) 0 ( 1 - P) = l + {d-l)P' and mean occupations at the low density sites of r r w - i P&R- = A g a i n the values of o here are probabilities of finding the left edge of a particle i n a particular lattice site. From these results we find the occlusion probability for the last site in the lattice is (1 + (d — I)(3)OR = (1 — (3). Finally, note that the currents and densities given in Eqs. 3.19 and 3.23 reduce to the known results for the standard T A S E P when d = 1. 3.2.4 The Matching Mean Field Phase Boundaries different approximations for the steady-state currents and densities cover all allowable values of the parameters a and 0. E q u a t i n g the current expressions given i n Eqs 3.14, 3.19, and 3.23, we find the following phase boundaries between the different current regimes ._, ^ . (Entry) a < a * 1 = — = — - T J Vd+l L a(l - a) —— l + {d-l)a = oL L = a l + (d-l)a a < f3 ( E x i t ) " * = VlTT < TTW^W °«= TW ^ 1 ^ J <- > 3 24 P <a 1 (Maximal) ot,B>—=—Vd+l 1 J m a x m a x = (—T=—— v^+l) 1 2 o~ = tMX max y/d(Vd+l) Note that for d = 1 the boundary and interior densities predicted by simple and refined mean field methods are equivalent. However for d > 1 we will see that the current and boundary densities in E q . 3.24 are significantly more accurate those predicted by the simplest mean field theory (Eq. 3.8). Eq. In fact, the currents and boundary densities of 3.24 match our M C simulation results to w i t h i n simulation error. T h e qualitative Chapter 3. T o t a l l y A s y m m e t r i c E x c l u s i o n Processes w i t h F i n i t e Size Particles 23 difference between the steady-state boundary densities predicted by E q . 3.8 and E q . 3.24 arise directly from the exclusion of the individual particles. Note that while the boundary densities are accurately predicted w i t h mean field theory, we expect that the predicted density profiles are as inaccurate as those produced by mean field theory for the standard d= 1 T A S E P [26]. A l t h o u g h the results i n E q . 3.24 were constructed by combining results from approx- imations unique to each of the three phases of the T A S E P , the predicted natures of the phase boundaries are consistent with the exact solution of the d = 1 T A S E P [15, 22, 26]. Specifically, ( d J / d a ) and ( d J / d f i ) are continuous i n the transition between the maximal current and boundary limited regimes, while the transition between the entry and exit limited regimes exhibits a discontinuity i n the first derivative of the current, and the boundary-limited to m a x i m a l current transitions have the following ^-dependent discontinuity in the second derivative (3.25) The d dependence of E q . 3.25 is different from that predicted using the simple mean field approach of E q . 3.9 and is shown below to compare favourably w i t h results from M C simulations. 3.3 Monte Carlo Simulations Extensive Monte Carlo simulations were performed to validate the various analytical models presented i n the previous two sections. A s we expected the densities in the lattice to vary significantly as the injection and extraction rates were varied, we based our Monte Carlo code on the B K L continuous-time algorithm [67]. T h e B K L algorithm was described in Chapter 2 and has the advantage of maintaining a constant computational efficiency over a wide range of particle densities. In the simulation, the particles themselves are set to occupy exactly one lattice site. However, all the particles were restricted to have a minimum inter-particle separation of d lattice sites. T h e m i n i m u m separation restriction Chapter 3. T o t a l l y A s y m m e t r i c E x c l u s i o n Processes w i t h F i n i t e Size Particles 24 was maintained i n the simulation by setting the rate of any particle movement that yielded an inter-particle separation of less than d sites to zero. E x c l u d i n g these modifications, the simulation proceeded exactly as described in Chapter 2 and A p p e n d i x A . Also, note that by modifying the B K L algorithm i n this way, the simulation exactly modeled the convention of defining the position of a size d particle by the location of its left edge. The magnitude of the finite size effect in our simulations was estimated by running with d = 1 and lattices of varying lengths. For lattices at least one thousand sites long, the M C results were found to systematically deviate from the known T A S E P results by less than half a percent. In the simulations w i t h larger particles, the length of the lattice was scaled w i t h d to ensure that at least one thousand particles could occupy the lattice simultaneously (N = lOOOd). T h e simulations were run for 4 x 10 steps, which was 9 sufficient to reproduce the known T A S E P results in lattices much longer than our d x 1000 site benchmark. In all simulations, the function R A N 2 i n [3], was used to generate the pseudorandom numbers i n the simulation. The period of this generator is approximately 2 x 10 3.3.1 18 and thus is much larger than the number of simulation events. Currents Figure 3.4 plots the steady-state current as a function of particle size d i n all three current regimes. W h i l e the simple mean field approximation always underestimates the currents, the refined mean field results given by Eqs. 3.19 and 3.23 give currents that match the Monte Carlo results w i t h i n simulation error ( « 1%). T h e currents for various values of d as functions of a and (3 are shown in F i g . 3.5. F r o m these curves one can estimate the critical values a* and (3* at which the phase boundaries lie. The mean field predictions for the currents displayed i n F i g . 3.5 are a combination of the results from two different equations. E q u a t i o n 3.14 was used for entrance and exit rates greater than the values of a*, (3* i n E q . 3.24. For smaller values of a, (3, Eqs. 3.19 and 3.23 were used. T h e simulated and refined mean field curves (Eq. 3.24) are symmetric w i t h respect to the interchange of a and (3. However, the simple mean field results (Eq. 3.8) always underestimate the currents and are shown for comparison by the Chapter 3. T o t a l l y A s y m m e t r i c E x c l u s i o n Processes 1 0.1 1 I I I I I I 0.075 t 0.05 w i t h F i n i t e Size • — -- Particles 25 M C Results Refined Mean Field Simple Mean Field (a) Entry Limited Current - - - _ 0.025 1 1 0.25 1 1 I I I I I I I I i i i i i i i i i i * * i i ~ 7 i ' (b) Maximal Current 0.2 - V 10.15 \ 0.1 '— -- _ 0.05 1 0.08 0.06 I I I I 4 5 I I • — ~* * I I i i \ - \ . \ i • -• • m i i i i ii (c) Exit Limited Current - \ 0.04 \ \ 0.02 0 1 ^ \ i 1 2 ~» i_ . . I 3 I I 6 I 7 i i i i i i 8 d 9 10 11 12 13 i 14 15 Figure 3.4: T h e currents w i t h i n each of the solution regimes as a function of particle size (d). Results from simple mean field, refined mean field, and Monte Carlo simulations are shown, (a) T h e entry limited current JL as a function of d for (3 = 10 and a = 0.1. (b) T h e maximal current J for a = (3 = 10. (c). T h e exit limited current J for a = 10 and (3 = 0.1. m a x R Chapter 3. Totally Asymmetric Exclusion Processes with Finite Size Particles 26 P Figure 3.5: Currents as a function of a or (3 for d = 2 ,4, and 8. (a) T h e entry limited to m a x i m a l current profile w i t h (3 = 10. (b) T h e exit limited to maximal current profile w i t h a = 10. T h e points are the Monte Carlo results, the solid lines are the refined mean field predictions (Eq. 3.24), and the dashed lines are the simple mean field predictions (Eq. 3.8). Chapter 3. T o t a l l y A s y m m e t r i c E x c l u s i o n Processes w i t h F i n i t e Size Particles 27 dashed curves. Note the clear asymmetry between entrance and exit rate limited currents in the simple mean field approximation. Figures 3.4 and 3.5 demonstrate the extremely good agreement between the M C results and the refined mean-field expressions (Eq. 3.24) for the current. T h e points in these plots each represent the mean of 10 M C runs of 4 x 1 0 events. T h e errors on the M C currents 9 were set equal to the standard deviation over the 10 M C runs, and were approximately one percent. In all cases, the difference between the optimal mean field predictions and the MC results was w i t h i n simulation error. T h e extremely high quality of the fit between the mean field currents and the M C results strongly supports our conjecture that the currents derived from the refined mean field approaches are exact. 3.3.2 Particle Densities We now compare the densities predicted i n Eqs. 3.8 and 3.24 w i t h those derived from MC simulations. T h e M C simulations also produce density profiles which are difficult to compute using mean field approximations. In Figure 3.6, the left boundary and interior densities, multiplied by d [phd and o d), max and the right boundary density a R are shown as functions of particle size d. T h e densities are the mean of 10 M C runs and have an error of approximately 3 percent. T h e densities or,d and cr , d, a n d o ma x R from our refined mean field analysis (3.24) match the Monte Carlo results exceptionally well. We have also verified that all the densities o^, c r m a x , and o R can be readily and ac- curately computed i n all current regimes. For example, oi can be found from aP{ji > d) = J, where J = J m a x for the maximal current regime, and J = JR for the exit limited regime. Similarly, the interior densities i n the entrance or exit limited regimes are found by equating E q . 3.13 to JL or J i n E q . 3.24 and solving for the appropriate root o. T h e R boundary and interior densities for all regimes are Chapter 3. T o t a l l y A s y m m e t r i c E x c l u s i o n Processes w i t h F i n i t e Size • o d (a = Q.l, t 1 1 2 3 1 4 5 1 1 6 d 1 7 1 8 28 (3= l ) L 1 Particles j 1 10 9 Figure 3.6: B o u n d a r y densities in each of the three current regimes. T h e solid lines show the predicted densities from the refined mean field approximations (Eq. 3.24). a l + (d- (Entry) (2d) -1 l + (d- l)a y/(l + (Exit) 1 d ad(l 0(1-0) + (d - 1)0) (2d)~ ^(i (Maximal) O-R C b u k l OL + a)/0 a(l - 1)J L l + (d- l)a (d-l)J ) -AdJ 2 L \ + L (d-l)J + l + (d-l)0 R (d-i)j y-4dj R R a (Vd+1) a d(^/d+l) 2 Vd(Vd+l) 0(Vd + l) ' (3.26) 2 where err, i n exit and m a x i m a l current regimes are the averages of the slowly varying densities over the first d lattice sites. Recall that i n the standard d = 1 T A S E P , mean field theory provides the correct boundary densities but incorrect density profiles [26]. Figure 3.7 depicts the occupation density profiles (o~d) as a function of site i for a T A S E P w i t h particle size d = 3. Notice that a weak anti-correlation arises at the right boundary even i n the entrance rate limited Chapter 3. Totally Asymmetric phase (inset of F i g . 3.7(a)). Exclusion Processes with Finite Size Particles 29 Similarly, there is a very small periodic behaviour in the density near the entrance in the exit rate limited regime (first inset of F i g . 3.7(c)). A more striking anti-correlation is evident near the right boundary of the exit limited lattice (second inset i n F i g . 3.7(c)) where the density profile splits into d branches with the branch corresponding to ON, aN-d, crN-2d, • • • having the highest occupation. T h i s is a consequence of the d-site exclusion strongly manifesting itself near the crowded exit, and our choice to represent the position of a particle by the location of its left edge. W h e n currents are exit rate limited, the last particle in the lattice spends a significant amount of time at its last possible site i = 1000<i T h e next site available for occupation is at i = 9999, thus forcing d—1 sites to be empty as long as the last particle occupies the last site. Correlations are small at the left boundary since a l l particles can only move to the right. A s particles are located by their left edges, and the observed correlations are due to close-packing of the particles near the lattice boundaries, i n partially asymmetric exclusion models we expect the inclusion of backward motion would strengthen the anti-correlation effect at the left boundary of the lattice, and weaken it at the right. 3.3.3 Phase Boundaries In order to map the phase diagram for the T A S E P w i t h finite sized particles, the boundaries between the exit limited, entry limited, and m a x i m a l current phases were determined. To find the phase boundaries from the Monte Carlo simulation results, current profiles in the a direction were computed for a range of fixed (3 values. Analogous profiles were computed in the exit limited phase. T h e following fifth order centred difference formula was then used to compute the first through third derivatives of the current data dJ dx (J(x 0 - 2 A x ) - 8 ( J ( x - A x ) - J(x 0 0 + Ax)) - J(x Q + 2Ax)). (3.27) Discontinuities i n the first derivative of J at the boundary of the exit and entry limited phases were marked by peaks in second derivative, while discontinuities in the second derivative of J at the boundary limited to m a x i m a l current transitions were marked by Chapter 3. T o t a l l y A s y m m e t r i c E x c l u s i o n Processes 0.5 w i t h F i n i t e Size Particles 30 (a) Entry Limited (a = 0.1, (3 = 1) 0.4 0.26 0.3 0.24 0.2 0.2^ 80 3000 (b) Maximal Current (a = 1, (3 = 1) 0.8 0.6 b~ 0.4 0.2 1.5 0.9 1 0.5 0.8 0 0 (c) Exit Limited (a = 1, (3 = 0.1) 1 2 0 500 10 2970 20 1000 1500 2000 3000 2500 3000 i Figure 3.7: Representative scaled density (aid) profiles i n the three current regimes, (a) T h e entry rate limited regime (a = 0.1,/? = 1). (b) T h e m a x i m a l current regime (a = B = 1). (c) T h e exit rate limited regime (a — l,@ = 0.1). T h e boundary regions are shown i n the insets. T h e strong splitting of the densities into high density and low density branches, is a product of the close-packing of the particles i n the exit limited regime, and the fact that particle positions are determined by the location of the particles' left edges. Chapter 3. Totally Asymmetric Exclusion Processes with Finite Size Particles 31 Figure 3.8: T h e phase diagram for the T A S E P w i t h d = 3 particles derived from the refined mean field results (Eq. 3.24) and the simple mean field results (Eq. 3.8). T h e points are the transitions estimated from M o n t e Carlo simulations. Chapter 3. T o t a l l y A s y m m e t r i c E x c l u s i o n Processes (a) w i t h F i n i t e Size Particles 32 M C Results '~~ Simple M e a n Field Refined Mean Field (b) 0.8 a = p (Refined Mean Field) a* (Simple Mean Field) *S (3* (Simple Mean Field) 0.4 0.2 5 d 10 6 Figure 3.9: (a) T h e discontinuity i n the second derivative of the steady-state current w i t h respect to a or 8 at the transition between the m a x i m a l current and the entry limited regimes. The simple mean field approach ( E q . 3.9) predicts an increasing discontinuity w i t h d, while the refined mean field approach (Eq. 3.25) predicts a j u m p that decreases w i t h d. (b) T h e critical values a* and 8* at which the boundary limited to m a x i m a l current transitions occur (Eqs. 3.8 and 3.24). peaks i n the t h i r d derivative. T h e transition point predicted by M o n t e Carlo was assumed to be at the value of a or 8 associated with the m a x i m u m value of these peaks. The phase boundaries predicted by the simple and refined mean field methods and Monte Carlo simulations can be viewed in F i g . 3.8. Additionally, and as can be seen i n Fig. 3.9, the M C results for the j u m p in the second derivatives d J / d a 2 2 and d J/df3 2 2 have noise magnified by the differentiation procedure. Nonetheless, the results agree quite well with the predictions from the refined mean field approaches. (Eq. 3.25). Chapter 3.4 3. Totally Asymmetric Exclusion Processes with Finite Size Particles 33 Summary We have presented mean field theories which accurately model the currents and densities of a totally asymmetric exclusion process w i t h particles of arbitrary size. Specifically, steady-state currents and densities were computed, using b o t h analytic, mean field approaches and Monte Carlo simulations. A l l calculations show that the steady-state currents qualitatively retain the three phase structure (low density, high density, m a x i m u m current) familiar from the standard (d = 1) T A S E P . Our best mean field model utilizes the canonical partition function for a set of size d particles confined to a finite lattice to produce current and density expressions in the maximal current regime. In the boundary limited regimes, an explicit enumeration of the occupation states near the lattice ends was employed. Unlike the simplest mean field approach, these refined analytic approximations provide simple formulae for the steady-state currents and densities for all particle sizes d that agree extremely well (within simulation error) w i t h our Monte Carlo simulation results. Based on this close agreement, we conjecture that the analytic expressions for the currents given by the refined mean field approaches (Eq. 3.24) are exact. B o t h analytic approaches predict first order transitions in the current between the entry and exit limited phases, and second order transitions in the current between the m a x i m a l current and the boundary limited phases. However only the refined mean field approaches produce estimates for the discontinuities in the current derivatives that match the simulation results. T h e extension of exclusion processes to larger particles (d > 1) allows one to directly apply the results to many biophysical problems including ribosome motion along m R N A [34, 69] and vesicle transport along microtubules. 34 Chapter 4 A Totally Asymmetric Exclusion Process with Periodic Structure 1 4.1 Background Continuing i n the theme of the preceding chapter, we now consider other modifications to the standard totally asymmetric exclusion process. A s described i n the introduction, the T A S E P successfully captures the discrete and driven nature of many physically interesting one-dimensional transport processes. However, in many systems the assumption of a single internal hopping rate p implicit i n the normal T A S E P does not adequately capture the full character of the transport behaviour. Examples include multi-step models for molecular motor motion, as well as models of ion transport through transmembrane channels w i t h selectivity filters and solvation zones. To improve the utility of the T A S E P i n modeling these sorts of systems, we generalize the open-boundary T A S E P to include two internal hopping rates. Previous studies of T A S E P s w i t h multiple movement rates focused on systems w i t h periodic boundary conditions [14, 31] and containing lattice sites where particles move with a reduced hopping rate. These studies show that the major effect of a single slow or "defect" lattice site was to induce phase segregation i n the periodic lattice with the formation of a region of rapid density change (a density shock) joining the two regions. In the case where multiple defect sites are randomly distributed, similar phase segregation can be seen and the current-density relation exhibits a noticeable plateau when the density is approximately 0.5. In open systems, the effect of a single slow or fast defect i n the interior of a lattice has been investigated [17]. W h e n the hopping rate associated w i t h the defect site was faster than the majority of sites i n the lattice there was no change T h e work in this chapter is based on "Steady-state properties of a totally asymmetric exclusion process with periodic structure", Greg Lakatos, T o m Chou, and Anatoly Kolomeisky, Phys. Rev. E 71 (2005) 011103 J Chapter 4. A T o t a l l y A s y m m e t r i c E x c l u s i o n Process 7" TT ? W t with Periodic Structure V 35 I Figure 4.1: T h e periodic dual-rate totally asymmetric exclusion model. A t every T lattice sites, the particle movement rate is p . A t all intervening sites, the movement rate is pi. A l l other aspects of the model are identical to those of the standard TASEP. 2 in the phase diagram, current, or bulk densities of the T A S E P compared to the uniform hopping rate case. Density profiles from Monte Carlo simulations do exhibit rapid density changes near the fast site, but these density perturbations rapidly decay and the density returns to the bulk value. In the case of slow defects, the phase diagram retains its three phase character, and the transition between the low and high density states remains on the a — P line. However, the transition between the boundary limited regimes and the maximal current regime become a function of the defect strength and approaches zero as the defect hopping rate is reduced. Additionally, in the m a x i m a l current state the particle current and bulk densities all become functions of the defect strength. Conversely, in the boundary limited regimes the currents and bulk densities are unchanged from the uniform hopping rate case. 4.2 M o d e l and Methods The generalization of the T A S E P investigated i n this chapter includes two internal hopping rates, p \ and p . Numbering the lattice sites starting from 1 at the far left, we assume 2 the p sites are arranged periodically with a period T and that the lattice consists of a 2 whole number of ( p i , p i , . •. ,p ) segments (Fig. 4.1). 2 We apply and compare two mean field approximations to the steady-state current and density of the dual-rate T A S E P in the maximal current, entry limited, and exit limited regimes. T h e accuracy of our mean field methods are determined through comparison with results from Monte Carlo simulations. Chapter 4.3 4. A Totally Asymmetric Exclusion Process with Periodic Structure 36 M e a n Field Theories 4.3.1 Simple Mean Field Methods We begin by considering the dual-rate T A S E P in the limit of large a and B, where the behaviour of the system will be determined entirely by the internal movement rates p\ and p . Ensuring continuity of the current i n a lattice w i t h periodicity T we find 2 J = p o (l 2 T - pi<7i(l <Ti) = - a) 2 = .. • = p i c r T _ i ( l - o) (4.1) T In writing E q . 4.1 we have assumed that the densities (oi) in the lattice have the same period T as the hopping rates i n the lattice. Solving E q . 4.1 for J i n terms of one of the densities o"j, then maximizing J, yields an approximation to the m a x i m a l current and densities in terms of p\ and p . Unfortunately, while it is easy to write E q . 4.1 2 for arbitrary T , the solution of E q . 4.1 yields increasingly unwieldy expressions for the current and the densities as T increases. A s a result we only show the expressions for the T = 2 case: j V\V2 = 2 ' *i = a = ^ T / - ' (4-2) /Pi 2 ^Pl + VP2 These relations show the expected invariance under the interchange of pi &ndp 2 Next we assume the dual-rate T A S E P retains the 3 phase character of the standard T A S E P , and look for approximations to the current and bulk densities in the boundary limited phases. To find the quantities i n the simple mean field approximation we assume a small value of a or B and assume the lattice contains a whole number of periods, so that the last lattice site is a p 2 site. We also assume that the steady-state densities in the lattice are periodic and have the same period as the hopping rates. Focusing initially on the entry-limited case and applying the current continuity conditions we find, for arbitrary T, Chapter 4. A T o t a l l y A s y m m e t r i c E x c l u s i o n P r o c e s s w i t h P e r i o d i c S t r u c t u r e J = a(l- <TI) = pi<ri(l - a ) = P20T{1 ~ •••= 2 37 (4.3) <7i). F r o m E q . 4.3, we find the following densities and currents for the T = 2 case Pi{p = _ - 2 a)a P2(Pi + o)-api ap (4.4) 2 p {pi + a) - api 2 o a>2 = Oil pi. Similarly i n the exit limited case we have J = P20"/v(l ~ °~N-T+l) = PlO~N-T+l(l - C/v-T+2) = . (4.5) P l C T j V l ( l - ON) = PO-N- G i v i n g , for T = 2 (p2-(3)i3 Pl P 2 ( p i +/?) - / 5 p i ' ^,1 1 - - , = (4-6) P2 = Pl(P2~ P) P2(Pl+P)~PlP' To determine the transitions between the m a x i m a l current and entry limited regions we equate the m a x i m a l current solutions generated by E q . 4.1 to the expressions for J a and J/3 in E q . 4.4 and E q . 4.6 and solve for the transition values a* and P*. For the T = 2 case this yields a =P a* = P2\JP\ r — VPI + — ( 4 . 7 , A ) VP^ F i n a l l y equating the current expressions i n E q . 4.4 and E q . 4.6 we find that the transition between the entry and exit limited regions occurs when a = 0. Chapter 4.3.2 4. A T o t a l l y A s y m m e t r i c E x c l u s i o n Process with Periodic Structure 38 Refined Mean Field Methods As we will see in the next section, the results of the simple mean field approach do not provide a particularly good match to the results of M o n t e C a r l o ( M C ) simulations. We attribute the poor performance of the simple mean field method to the method's failure to capture the correlations between the occupation probabilities of different lattice sites. To partially address this deficiency, we develop a more refined mean field approach. T h e basis of the refined mean field approach is to assume that the density in the T A S E P lattice exhibits the same periodicity as the hopping rates. We then allow correlations between lattice sites w i t h i n a single period, and neglect correlations between adjacent periods. T h i s approach is related to the cluster mean field methods developed for other non-equilibrium stochastic lattice gas models [70-72]. We begin by applying variations of this mean field approach to the T = 2 case where we will be able to produce concise and accurate expressions for the currents and densities i n all three phases of the dual-rate TASEP. We first consider the m a x i m a l current phase of the T = 2 dual-rate T A S E P , and begin by defining the pair probability P ( x i , x i ) . i + T h e pair probability P ( x j , x i + 1 ) is the probability of finding a p i site w i t h occupancy X ; , followed by a p site w i t h occupancy 2 x i. i+ T h e time evolution of the occupancy state of any two adjacent sites i n a T A S E P will depend on the two sites themselves along w i t h the pair of lattice sites immediately to the left or the right of the two site grouping. Thus the master equation for the two site probability P ( 0 , 0) is ' K = dt -p (P(0,l,0,0) + P(l,l,0,0))+p (P(0,l,0,0) + P(0,l,0,l)). 2 2 (4.8) Now assume that each pair of ( p i , p ) sites behaves like a statistically independent unit, 2 uncorrelated w i t h the other (p\,p ) 2 the probabilities i n E q . P ( X J , x i, i+ x , i+2 x) i+3 Using these assumptions 4.8 can be decomposed into products of pair probabilities, « P(xi, x P ( 0 , 0) in the steady-state: periods in the lattice. i + i)P(x i + 2 , x ). i+3 T h i s yields the following equation for Chapter 4. A Totally Asymmetric = Exclusion Process with Periodic p (P(0,l) -P(0,0)P(l,l)) = 0. 2 2 Additionally, given the definition of P(xi,Xi+i), = P(1,0) + P(1,1), 0 = P(0,1) + P(1,1), P(0,0) = 1-P(0,1)-P(1,0)-P(1,1). 2 39 ( 4 9 ) the following relations hold 0 1 Structure (4.10) Using Eq. 4.10 in Eq. 4.9, produces the relation 7^'. 1 + ( c r - O-i) (4.11) 2 An equation for (5(1,0), the probability of having an occupied p site followed by an 2 unoccupied pi site, can be found from Eq. 4.11 by interchanging p and p , as well as o\ 1 and o . 2 Applying the current continuity condition P i P ( l , 0 ) = p Q(l,0), 2 2 results in the relation o-i(1- cry) - (02 ~ oif 1+ (fJ 2 - _ o {l - Qi) - ( " i - 02? 2 1 + <7l) (<Ti - <7 ) 2 Now consider the substitution o — 0\ = /c and assume pi to be the larger of the two rates. 2 Thus oi < o and k > 0. Solving Eq. 4.12 for o\ yields 2 o l = (2k( Pl + p ) - 2{ 2 - Pl p )Y l 2 ( p ( l - k ) - (k 2 2 Pl - l) ) + 2 (4-13) V ( l - k ){k( +p ) 2 Pl 2 - ( i - 2)){Sk( P P Pl + p) 2 Equation 4.13 shows that o\ is real when — 1 < k < 3 ^ ^ ) or Eq. 4.13 into Eq. 4.12 gives J = -, 1 1 fc> for 0 < k < 2i-jE2, and § P1+P2' dk fJ k lV Pl 2 < k < 1. Substituting , >, which shows J < 0 for k > (Pl-P2)-«(P1+P2) ' = { -p )) P I P ^ - P ' ) J > 0 P1+P2' > 0 V k. As J must be positive, and (fc(p2+pi)+(P2-pi)r <7i must be real, and J is monotonically increasing with k, the maximum value of the current must occur when k = 3^"+^) • Substituting this value of k into Eq. 4.13, yields Chapter 4. A Totally Asymmetric Exclusion Process with Periodic Structure 40 the following expressions for the current and bulk densities in the maximal current phase Pi + 2p 2 3(pi + p ) ' 2 2 p i + P2 0"2 3(pi+P2)' J ( - ) 4 14 P1P2 = 2(pi+p )' 2 Note that the solution is invariant under the exchange ( c r ^ p i ) for ( c r , p ) , and reproduces 2 2 the standard T A S E P results in the limit pi = p : 2 We can use the results of Eq. 4.14 to predict the location of the phase boundary between the maximal current and entry limited phases. Using the results for the densities in the maximal current phase and applying current continuity at the entrance of the lattice we find " M 1 l -S( ^+p ^)J l Pl = n ^ -)v 2(^+p 2 Pl (4-15) v ; 2 Giving a value of a = 2 p ^ y ( 4 - 1 6 ) To find the critical value of the extraction rate (3*, we enforce current continuity at the exit of the lattice, yielding F = P 3PlP2 2(2pi+p )- (4 17) 2 An analogous approach can be used to find current and density estimates in the entry and exit limited regions. Considering the entry limited case first, we find the following master equation for the occupancy of the first two sites in the lattice: Chapter 4. A Totally Asymmetric rfP(0,0) dt riP(0,l) Jt Exclusion Process with Periodic Structure = - a P ( 0 , 0) + p P(0,1)(P(0,0) + P(0,1)), = -c*P(0,1) - p P ( 0 , 1 ) ( P ( 0 , 0 ) + P(0,1)) 41 2 2 +PiP(l,0), (4.18) dP(l,0) dt = aP(0,0) + p P ( l , 1)(P(0,0) + P(0,1)) 2 -PiP(l,0), dP(l, 1) dt = - p P ( l , 1)(P(0,0) + P(0,1)) + « P ( 0 , 1 ) . 2 Similarly, in the exit limited phase we find the following master equation for the last two sites in the lattice rfP(0,0) /5P(0,l)-p P(0,0)(P(0,l) 2 dt dP(0,l) + P(l,l)), -/3P(0,1) - p P(0,1)(P(0,1) + P ( l , 1)) 2 dt +PiP(l,0), rfP(l,0) ( . ) 4 19 /?P(1,1)+ P(0,0)(P(0,1) + P(1,1)) dt P2 -PiP(l,0), dP(l,l) dt - -/3P(1,1) +p P(0,1)(P(0,1) + P ( l , 1)). 2 Equations 4.18 and 4.19 assume that the occupancy probabilities P(xi,x i) i+ do not vary significantly as a function of position near the entrance and exit of the lattice. This assumption is analogous to that used with the standard T A S E P in the boundary limited regions, where the density is assumed to be essentially constant near the rate limiting boundary, and is supported by the results of Monte Carlo simulations for the T = 2 T A S E P (c. P(xi,x i) i+ f. Fig. 4.6). Applying the normalization condition on the probabilities and solving Eqs. 4.18 in the steady-state gives two sets of solutions for the Chapter 4. A T o t a l l y A s y m m e t r i c E x c l u s i o n Process with Periodic Structure 42 current and bulk densities in the entry limited phase. We choose the solution which gives o~\ —• 0, c —> 0 in the limit a —> 0. This yields the following results for the current and 2 bulk densities in the entry limited phase « \P1P2 - a ( p i + a) + y/4pla(p2 - a) + (pi(p - <*) + a ) ) 2 2 2 2p (Pi + a) 2 a p + (a + pi)(a + p ) - y / b p j a i p i _ a ' 2 2 - a) + (pi(P2 - a) + a ) 2 (4.20) 2 2p (p! + a) X 2 a — • = Ca,2 P2 Similarly, solving Eqs. 4.19 and choosing the solution that gives u\ —> 1, cr —> 1 in the 2 limit P —+ 0, produces the following approximations to the current and bulk densities in the exit-limited phase P [PlP2 - P(Pl +P)+ J P V4P(P(P2 + (Pl(P2 ~P)+ P) 2p ( +P) = 2 *M ~P) = 2 2 ' PL 1 - A ( - ) 4 P2 PXP2 ~ P(Pl +P) + VMPJP2 ^ ~P) + (Pl(P2 ~ 21 P ) + W 2p ( +/3) 2 Pl Examining Eqs. 4.20 and 4.21, it can be seen that the refined mean field method predicts that the current is invariant under the exchange of a and P, and that (1 — cr ) = o p ^ , a>1 and (1 — <7 a)2 ) = o/3 i t under the exchange of a and p . These relationships are a product of the underlying particle-hole symmetry in the T = 2 TASEP, and will be discussed in more detail later in the chapter. In principle the method used to generate Eqs. 4.20, and 4.21 could be used to generate analytic expressions for the currents and densities in all three phases of the dual-rate T A S E P with an arbitrary period T. However for T > 2, the results rapidly become unwieldy. As a result, we seek an automated means of producing these estimates with the expectation that the mean field results will be reasonably accurate and will require less Chapter 4. A T o t a l l y A s y m m e t r i c E x c l u s i o n Process 2 43 •'2 P 4 _ with Periodic Structure °eff P,eff 1=1 1=2 r=3 1 = 4 L = 4 Segment Figure 4.2: The arrangement of lattice sites used for an L = 4 finite segment mean field (FSMF) approximation to the T = 2 TASEP. The master equation for the state of the four marked lattice sites is solved exactly. The four sites are coupled to the rest of the lattice by assuming an effective injection rate of cieff = P2O4 on the left, and an effective extraction rate of /3ff = p ( l ° i ) the right. — e o n 2 computer time than running a Monte Carlo simulation. To begin, we focus on developing the method for computing the maximal current phase properties. Consider a finite lattice containing L sites and a whole number of ( p i , p i , . . . , p ) 2 periods (Fig 4.2). Assuming the bulk densities in the T A S E P lattice are periodic and inherit the same periodicity as the movement rates, we define an effective injection rate (ct ff) and an effective extraction rate (/? ~) for this finite segment. Specifically, we set e e "eff = P2OL and /3ff = p ( l — (7i). Using these rates and the known values of p\ and ~>, e 2 2 we can build the transition matrix for the master equation describing the motion of the particles in the L-site finite segment. Recall that the Master equation can be written dP{t) dt MP(t). (4.22) Here M is the matrix of transition rates that describe the T A S E P and P(t) is a vector of the time dependent probabilities of finding the T A S E P lattice in a particular occupation state. The transition matrix M is typically sparse, and even for a small segment it is also large (2 L x 2 ). To construct M efficiently it is desirable to avoid having to test the L connectivity of every possible state of the finite segment. This can be accomplished by first defining a vector r containing the internal hopping rates of the finite segment, and supplying an initial guess for the densities o\ and or,. Specifically, set the i equal to the hopping rate of the i t h t h entry of f lattice site in the finite segment. The hopping rate of the last site in the finite segment however is set to /?fr> so that r e L — B . With the finite e f i segment and the rate vector (f) defined, the algorithm to construct the master equation Chapter 4. A Totally Asymmetric Exclusion Process with Periodic Structure 44 Figure 4.3: T h e steps i n the algorithm to generate the transition matrix for a T A S E P . For the purposes of illustration a three site lattice has been used. Each possible occupancy of the lattice is treated as a bit pattern, and each state is labelled w i t h the corresponding decimal value (i.e. lattice occupancy O i l —> state 3). T h e states are divided into two groups; states where the first lattice site is occupied (1-states), and states where the first lattice site is empty (0states). Regardless of the number of lattice sites i n the T A S E P , the transitions between the two classes of states always occur between the first half of the 1-states where the second lattice site is unoccupied and the second half of the 0-states where the second lattice site is occupied, (solid arrows i n the figure). Calling the algorithm recursively on both the 1-states and the 0-states, and ignoring the highest order bit (i.e. the leftmost lattice site), enumerates all the remaining state transitions (dashed arrows). F i n a l l y the transitions between each 0-state and each 1-state generated by injection at the left edge of the lattice (dotted arrows) are enumerated. transition m a t r i x is (see F i g . 4.3) 1. L a b e l t h e o c c u p a n c y s t a t e s o f t h e Label each occupancy state of the finite segment finite-segment w i t h a number determined by treating the occupancy of the lattice as an integer value expressed in base-2. For example, an L = 3 lattice w i t h only the last site filled (i.e. an occupancy of {001}), is said to be i n state 1, while the same lattice w i t h the last two sites filled (an occupancy of { O i l } ) is said to be in state 3. 2. S e p a r a t e t h e s t a t e s i n t o t w o g r o u p s Separate the states labelled i n step 1 into two groups; one group containing the Chapter 4. A Totally Asymmetric Exclusion Process with Periodic Structure 45 states with the left-most lattice site unoccupied (the 0 - s t a t e s ) , and the other group containing the states with the left-most lattice site occupied (the 1 - s t a t e s ) . 3. E n u m e r a t e the internal transitions between the 1-states and the O-states Excluding particle injection, the transitions between the 1-states and the 0-states will always occur between the lowest half of the 1-states, and the highest half of the 0-states. For example, in an L = 3 finite segment, the transitions take state 4 ({100}) to state 2 ({010}), and state 5 ({101}) to state 3 ({011}). To determine the rate of the transitions, we make use of the rate vector r. Specifically, we consider the current recursive iteration of the algorithm, and use this as an index for the rate vector. For example, during the second recursive iteration of the algorithm, all the enumerated transitions will involve the movement of a particle from site 2 to site 3. Thus, all the transitions occur at rate r ; the rate associated with the second 2 site in the finite segment. In general, the i t h recursive iteration of the algorithm will enumerate transitions where a particle occupying site i of the finite segment moves to site i + 1. These transitions occur at rate r*j. 4. Recurse on the 0-states and on the 1-states Recurse on the 0-states, and the 1-states, ignoring the left-most lattice site and starting the recursion from step 2. For example, if on iteration i we are working with an L = 3 site segment, in iteration i +1, we work with two L = 2 site segments. The first of these L = 2 site segments is generated by ignoring the left-most lattice site of the 0-states of the L = 3 segment, while the second L = 2 site segment is generated by ignoring the left-most lattice site of the 1-states. 5. Enumerate the injection events Once the recursion specified in steps 1 through 4 has completed, enumerate all the transitions representing injections into the finite segment. These transitions take the 0-states to the corresponding 1-states at rate a rf. For example, if we are working e on an L = 3 finite segment, we enumerate the following transitions; 0 ({000}) —» 4 ({100}), 1 ({001}) - 5 ({101}), 2 ({010}) - » 6 ({110}), and 3 ({011}) - 7 ({111}). Chapter 4. A Totally Asymmetric Exclusion Process with Periodic Structure 46 All these transitions occur at rate a ff. e Using this recursive algorithm to produce the transition matrix for the finite-segment master equation, we can apply an iterative procedure to find the densities and currents through any finite segment of the T A S E P lattice. This iterative procedure proceeds as follows 1. B u i l d the M a s t e r E q u a t i o n Given an initial guess for the values of (<TI,CTL),use the recursive algorithm to produce the transition matrix (M(<7i, OLIVIIVI)) The finite segment densities (OI,OL) a r e f ° the finite-segment master equation. r introduced into M solely through the ef- fective injection rate ct s = P2&L, and the effective extraction rate /? ff = p ( l ~~ e e 2 2. Solve the M a s t e r E q u a t i o n To find the steady-state currents and densities, compute the eigenvector (V(ci,o , i 2 P )P2 of M(CTI, c r , p , p ) associated with the zero eigenvalue i 1 2 M(o ,o , ,p )V 1 2 Pl = 0. 2 (4.23) Normalize V using the expression (4.24) i The vector element Vi gives the steady state probability of finding the finite segment in the state with label i. For example, with an L = 3 finite segment, V3 gives the steady state probability of finding the finite-segment in occupancy state {Oil}. 3. C o m p u t e new Density Estimates With the steady state occupancy probabilities from step 2, compute new estimates for U\ and or, Chapter 4. A Totally Asymmetric Exclusion Process with Periodic E Structure 47 > v (4.25) 4. Start a new iteration or finish Compare the new density estimates (ol,o* ), to the previous estimates L If (OI,OL)- |cr* — o !| < e and \o* — OL\ < e, then the iteration has reached a fixed point and - L the procedure ends. Here e is an arbitrary convergence parameter (set to 10~ in 4 this study). If the convergence criteria are not met, return step 1, setting o\ = cr*, and OL = o* . Note that as only the effective injection and extraction rates change, L a complete reconstruction of M is not required. Once a fixed point in the finite-segment densities (O~\,O~L) is finally reached, the current through the finite segment can easily be computed from the expression J = a ( l - o-i) = P2<y {l eff L oi). (4.26) For convenience we refer to the use of these two algorithms to compute the steady-state properties of a T A S E P as the finite segment mean field (FSMF) method. Thus far we have described an implementation of the finite-segment method useful for generating current and density estimates in the maximal current phase. However, the method can be extended to treat the boundary limited phases as well. This can be accomplished by setting Qf rf = a or (3 $ = (3 as appropriate, and applying the density e e self-consistency condition at the end of the finite-segment that lies in the interior of the T A S E P lattice. While a simple and effective approach that can be generalized for use with variety of defect arrangements in addition to the periodic arrangements explored here, the exponential increase in the size of the transition matrix with the length of the segment quickly renders Eq. 4.23 analytically intractable. As a result, Eq. 4.23 is solved numerically. Employing the eigensolvers in the well-known A R P A C K linear-algebra software library [1] and built into the commercial program M A T L A B , we could treat segments containing Chapter 4. A Totally Asymmetric Exclusion Process with Periodic Structure 48 < 18 lattice sites. 4.4 Monte Carlo Simulations Monte Carlo simulations were performed to validate the analytic and numerical models presented in the previous section. As we expected the densities in the lattice to vary significantly as the internal hopping rates were varied, we based our Monte Carlo code on the B K L continuous-time algorithm [67]. As was the case in Chapter 3, in the dual-rate TASEP, the B K L algorithm has the advantage of maintaining a constant computational efficiency over a wide range of particle densities. To simulate the dual-rate TASEP, the only modification made to the algorithm outlined in Chapter 2 was to set the rate for all transitions involving the movement of a particle located at a lattice site with an index i — s T + 5 equal to p 2 The magnitude of the finite size effect in our simulations was estimated by running lattices of varying lengths. For lattices one thousand sites long, the M C results were found to systematically deviate from the known T A S E P results ( = p ) by less than half a Pl 2 percent. As a result, unless otherwise noted, we used lattices containing 1000 periods for all our simulations. The simulations were run for 7 x 10 steps, which was sufficient 8 to reproduce the known T A S E P results in lattices much longer than our 1000 period standard. In all simulations, p \ was normalized to 1. 4.4.1 Currents The maximal current results (Fig. 4.4(a)) show the expected qualitative behaviour with a low current for small values of p 2 and a value of 1/4 when p 2 = Pi = I- Note that each M C result displayed in Fig. 4.4 is the mean over 10 separate runs. The error was assumed to equal the standard deviation of the set of 10 current results and is approximately 1% of the displayed values. Consistent with our expectations, we find that the refined mean field approach, and the finite-segment mean field approach, provide better approximations to the Monte Carlo results than does the simplest mean-field method. The relatively poor Chapter 4. A Totally Asymmetric (a) 0.25 0.2 • — --• — Exclusion Process with Periodic Structure 49 M C Results Simple mean field Refined Mean Field L = 12 FSMF ^ 0.15 0.1 0.05 0 (b) 0.25 0.2 ^ 0.15 0.1 0.05 Q 0.4 0.6 "2 Figure 4.4: (a) T h e Monte Carlo ( M C ) , simple mean field, refined mean field, and L = 12 finite-segment mean field ( F S M F ) predictions for the T = 2 T A S E P current in the m a x i m a l current phase (a = 10, (5 = 10). (b) M a x i m a l currents for dual-rate T A S E P s of various periods. T h e F S M F estimates shown by solid lines were produced using a single period of the lattice as the finite segment. W h i l e providing reasonable estimates, the deviation between the M C (a = 10, 6 = 10) and F S M F (lines) results increases as the period (T) increases. However (c) displays the percent difference (() between the Monte Carlo and F S M F results for a T = 5 T A S E P , and shows that the increase in the F S M F error w i t h increasing T can be partially reduced by increasing the segment length. Chapter 4. A Totally Asymmetric Exclusion Process with Periodic Structure 50 performance of the simple mean field model can be ascribed to the pair correlations present in the two-rate TASEP, which the simple mean field model completely fails to capture. The refined maximal current mean field results (Eq. 4.14) clearly performs best when p/i 2 P ~ 1. In this limit the model is essentially the normal single rate TASEP, and the expected pair correlations between adjacent p\ and p sites are small. 2 Figure 4.4(b) shows the maximal current predicted by both Monte Carlo simulations and the F S M F method in TASEPs with different periods, as a function of the slow hopping rate p . The F S M F method clearly performs better for smaller values of T. For larger T , 2 the method underestimates the Monte Carlo currents for 0.1 < p < 1.0 and overestimates 2 for very small values of p . 2 As shown in Fig. 4.4(c) the error in the F S M F predictions can be reduced by increasing the number of periods in the finite segment. However the improvement in the finite segment estimates grows sub-linearly, while the computation cost of the method increases exponentially. As the F S M F method directly solves the Master equation for the L site finite segment, the error in J is partially due to the meanfield coupling of the finite segment to the remainder of the lattice. The fact that the finite segment underestimates the Monte Carlo current for p > 0.1 suggests the existence 2 of an anti-correlation between the p site of one period and the first pi site of the next. 2 Additionally, the fact that the F S M F estimates improve with increasing L also indicates that there may be long-range correlations that span multiple periods and thus cannot be accurately captured when L = T. The current predictions for the entry-limited regime are shown in Figs. 4.5(a) and (b). While not shown, the current predictions for the exit-limited regime are identical, and the reason for this will be described later in the chapter when the particle-hole symmetry of the two-rate T A S E P is addressed. Figure 4.5 shows that the F S M F accurately captures the growth in the current as a is increased over a wide range of T values. However deviations between the M C and F S M F results increase as the transition to the maximal current regime is approached. Note that when maximized with respect to a or /?, the boundary limited variants of the F S M F method produce currents that match those produced by the maximal current variant. Chapter 4. A Totally Asymmetric 0 Exclusion Process with Periodic Structure 51 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 a Figure 4.5: (a) Current profiles along the a direction in the T = 2 dual-rate T A S E P phase plane (P = 10). The points are Monte Carlo results, the dashed lines were produced using Eq. 4.4, and the solid lines were produced using an N — 8 F S M F approach. Eq. 4.20 and the F S M F approximation are indistinguishable, (b) Current profiles for dual-rate TASEPs of various periods with p = 0.25. The F S M F results (lines) were produced using an L = T F S M F approach (/? = 10). 2 Chapter 4. A Totally Asymmetric Exclusion Process with Periodic Structure 0.6 W) 0.5 0.4 0.3 0.2 0.1 Entry Limited (a - 0.05, (3 = 10) (b) Maximal Current (a = 10, (3 = 10) 0.9 [(c) 0.8 0.7^ 0.6 \ 0.5 P 0.4 0 Exit Limited (a = 10, (3 = 0.05) 0.8 52 .0.6 0.4 0.2 500 1000 i 1500 2000 Figure 4.6: Monte Carlo density profiles for a T — 2 T A S E P with p = 0.1. Comparison with Fig. 2.1 shows that the profiles in the T = 2 T A S E P retain the same qualitative character as the profiles in the standard T A S E P . Note that in the boundary limited regimes the density profiles are approximately constant near the rate limiting boundary, consistent with the assumptions leading to Eqs. 4.18 and 4.19. 2 4.4.2 Densities Referring to Figs. 4.6 and 4.7, we find the finite-segment mean field, and the refined mean field methods, provide excellent matches to the densities produced by Monte Carlo simulations. The densities produced by the refined mean field method, and the finite segment method are all within 5% of the M C results. The M C results themselves represent an average over i P and p 2 sites in the centre of the T A S E P lattice. The mean density at the two classes of sites in the T A S E P lattice is then averaged over the results of 10 separate runs. The error on the displayed densities is approximately 3%. For both mean field approaches, the quality of the agreement is relatively uniform over all the values of Figure 4.8 displays the density profile in the centre of the lattice for TASEPs with T = 5 and T = 9. Far from the boundaries, the density profiles show the expected Chapter 4. A Totally Asymmetric 1 i ii iii Exclusion Process with Periodic Structure — T "— 1 — 1—I — 1 0.8 I I I 53 I Refined Max Current MF - - Simple MF - - L = 4 FSMF ,0.6 ef 0.4 0.2 P/ / / i Q 0 0.2 0.4 0.6 Pi _l 1 I II I ..I. 1 . . 1. 0.8 i. I_t- Figure 4.7: Maximal current phase (a = 10,/? = 10) densities for a T = 2 T A S E P from Monte Carlo simulations, the refined maximal current mean field method (Eq. 4.14), the simple mean field method (Eq. 4.2), and an L = 4 finite segment mean field approach. The Monte Carlo results are averages over the central 50 sites of a 2000 site lattice. The simple mean field assumption produces the largest errors at small values of p where the density correlation between adjacent sites is expected to be large. 2 periodicity in all three phases. While the F S M F results show good agreement with the M C results, there is a tendency to slightly underestimate the densities in the exit limited phase, and slightly overestimate in the entry limited phase when p = 0.2. In the maximal 2 current phase the profiles in any one period are approximately linear which is similar to the behaviour of a T A S E P lattice on the boundary between the entry and exit-rate limited regimes. Finally, Fig. 4.9 shows the densities as a function of a in the entry-limited regime. While the F S M F results track the growth of the densities with increasing a there is a tendency to overestimate the densities when p ~ 0.5. The exit-limited results are related 2 to the results displayed in Fig. 4.9 by particle-hole symmetry which will be discussed later in the chapter. Chapter 4, A Totally Asymmetric -10 -5 Exclusion Process with Periodic 0 5 Structure 54 10 i Figure 4.8: The bulk density profiles for the T = 5 and T = 9 dual-rate TASEPs with p = 0.2, in the maximal current (a = 10,/? = 10), entry rate limited (a = 0.05, P = 10), and exit-rate limited (a = 10,/? = 0.05) phases. The plots show the good agreement between the Monte-Carlo (points) and L = T F S M F (lines) results. The M C results are the central 5 periods (T = 5), and central 3 periods (T = 9) of the lattice. 2 Chapter 4. A Totally Asymmetric Exclusion Process with Periodic Structure 55 a Figure 4.9: The density profile along the a direction in the T = 2 dual-rate T A S E P phase plane ((3 = 10). The dashed lines were produced using Eq. 4.4, while the solid lines were produced by an L = 8 finite-segment mean field approach. The Monte Carlo results are averages over the central 50 sites of the T A S E P lattice. The F S M F results are indistinguishable from those of Eq. 4.20. Although the solution for o is the same in Eqs. 4.4 and 4.20, the predicted value of a* differs between the two mean field theories. This is the reason for the significant difference between the o profiles predicted by the two mean field methods at small values of p . 2 2 2 Chapter 4.4.3 4. A Totally Asymmetric Exclusion Process with Periodic Structure 56 Phase Diagram The phase diagram derived from simulations is displayed in Fig 4.10(a) for the T = 2 case, along with the phase boundaries predicted by the various mean field methods. The Monte Carlo phase boundaries were computed by taking numerical derivatives with respect to a and ft of the mean value of a in the centre of the lattice, and locating any clear discontinuities in the derivatives. The dual-rate T A S E P retains the general 2 form of the standard T A S E P phase-diagram, and the order of the transitions remains unchanged; first order in the current between the high and low density regimes, and second order between the boundary limited and maximal current regions. The most significant deviation from the standard phase diagram is the increase in the maximal current region which accompanies a decrease in one of the hopping rates (p in our example). Physically, 2 the maximal current region is defined as the region of the (a, (3) parameter space where the internal motions determine the net particle flow through the lattice. The increase in the area of the maximal current phase with decreasing p is then expected, as decreasing 2 p slows the motion of the particles in the interior of the lattice. 2 Despite the varying degrees of success in predicting accurate steady state currents and densities, all three mean field approaches predict the phase boundaries within approximately 20% error. Note that the phase transitions predicted by the boundary-limited results were determined numerically by maximizing J in Eqs. 4.20 and 4.21. Figure 4.10(b) displays the change in the phase diagram as the period T is increased and p 2 = 0.1 is held constant. The major effect of increasing T at a fixed p is to increase the values of 2 a* = ft*. 4.5 Other Periodic Arrangements The dual-rate T A S E P model investigated in the previous sections was built on lattices composed of an integer number of {p\,Pi, i, P • • • ,p } periods. While a useful extension to 2 the TASEP, it is also interesting to investigate other periodic arrangements of defect sites. For example we can consider the case where the first defect is located at an arbitrary site 5 Chapter 4. A T o t a l l y A s y m m e t r i c E x c l u s i o n Process with Periodic Structure 57 Figure 4.10: (a) The phase diagram for the T = 2 two-rate TASEP. The dashed lines are produced by the simple mean field theory (Eq. 4.7). The dotted lines are produced by the maximal current results (Eqs. 4.16, 4.17), and the solid lines are produced by the boundary-limited results(Eqs. 4.20, 4.21), while the points are the Monte Carlo results, (b) The phase diagram for the tworate T A S E P at various periodicities, with p = 0.1. To simulation accuracy, the transition between the high and low density phases occurred along the a = (3 line for all the periodicities displayed. Points are Monte Carlo results, lines are drawn as a guide to the eye. 2 Chapter 4. A Totally Asymmetric Exclusion Process with Periodic Structure 58 (a; 6 = 1 a '1 9 (b) 8 = 2 a 8 (c) Pi 5 = 3 -T—- a I] T Figure 4.11: Variations on the periodic dual-rate TASEP. The three lattices displayed correspond to (a) 5 = 1, (b) 5 = 2, (c) <5 = 3. In each case the lattice length was chosen to ensure the resulting T A S E P possessed particle-hole symmetry. within the first T sites of the lattice, as well as lattices that do not contain a whole number of periods (Fig. 4.11). When considering these variations on the dual-rate T A S E P model however, the issue of particle-hole symmetry arises. To investigate the conditions under which a periodic dual-rate T A S E P would preserve particle-hole symmetry, first define a vector f consisting of all the movement rates in the T A S E P lattice. Specifically, r is constructed so that element r, gives the rate at which particles in lattice site i move to the right. Additionally, N is defined to be the total length of the T A S E P lattice. Now consider the probability P(a, 6, r, xi, x 2) the T A S E P lattice in an occupancy state xi,x , 2 • • •, XN) of finding • • • ,x.jv with injection rate a, extraction rate /?, and hopping rates r. Here x. = 1 if lattice site i is occupied, Xi = 0 otherwise. h Particle-hole symmetry then requires P(a,f3,f,{xi,x ,...,x }) 2 where = P(B,a,f,{x ,x -i,...,x }) N N N (4.27) l = 0 if Xj = 1, and x — 1 if x = 0. For Eq. 4.27 to hold, the rate of the particle x x movement that takes the T A S E P lattice from a occupancy state s\ = { x i , x , . . . 2 to an occupancy state s must equal the rate of the particle movement that takes the lattice 2 from state S\ = {x ,x _i,... N N ,x{\ to occupancy state s . This places constraints on the 2 possible values of A f that can produce particle-hole symmetry. To determine the symmetry-preserving values of N, consider a lattice occupancy state Chapter 4. A Totally Asymmetric Exclusion Process with Periodic Structure 59 S i where a particle located at site k moves to site k + 1, changing the lattice state to s - By definition this movement would occur at rate r^. In lattice state si the equivalent 2 move involves a particle hopping from site N — k to site N — (k — 1) at a rate rjv_fc. The definition of the particle-hole symmetry operation (Eq. 4.27) ensures that the final state is then s~2. To have these exchange-symmetric moves occur at the same rate, equal r u - h - must = p . Then by definition of the dual rate Now consider the case where 2 TASEP, k = m T + 5 for some integer m (see Fig. 4.11). Particle hole-symmetry will then be satisfied if N - m T - 8 = n T+ S (4.28) => N [n + m ) T + 28 for an arbitrary integer n. Thus the dual-rate T A S E P model will satisfy particle-hole symmetry when N = w T + 28 for some integer w > 0. Geometrically, satisfying Eq. 4.28 guarantees that if the first defect is located at site 8 of the lattice, the last defect will be located at lattice site N — 8. To see that a value of N satisfying Eq. 4.28 will never confuse particle movements occurring at rate p \ for those occurring at rate p , consider the case when k = m T + 7 2 where 7 < T , 7 7^ 8. Under symmetry exchange, a particle move occurring at rate would then be mapped to a move occurring at rate TN-k = (4.29) r(u!-m)T+6+(5-y) As 7, 8 < T and 7 7^ 8, there are no integers m , n, w such that (w — m ) T + 8 + (8 — 7) = n T + 8 . Thus a particle movement occurring at a rate p i cannot be mapped to a movement occurring at rate p under the symmetry exchange operation. A n interesting consequence 2 of the restriction on is that a lattice with the first p site located at site T and containing 2 a whole number of periods preserves symmetry regardless of the value of T, while a lattice with the first p site at lattice site 1 and containing a whole number of periods violates 2 symmetry if T ^ 2. Figure 4.12 shows the phase diagrams for various symmetric and non-symmetric realizations of the dual-rate T A S E P with T = 3 and p = 0.25. The results displayed in 2 Chapter 4. A T o t a l l y A s y m m e t r i c E x c l u s i o n Process with Periodic Structure 1 (a) - 5 - 5 - 5 0.8 0.6 ft 1,N = = 2,N = = 3,N = = \,N = = 2,N = = 3,N = = • 5 • 5 • 5 60 1205 1204 1203 1203 1203 1204 CO. f 4 0.4 0.2 0.2 0.4 0.6 0.8 cx Figure 4.12: Phase diagrams for the T = 3 dual-rate T A S E P for various values of 8 and N derived from Monte Carlo simulations (p — 0.25). The lines are the best-fit phase boundaries for TASEPs with values of that satisfy Eq. 4.28. The points are the Monte Carlo predictions for the phase boundaries of TASEPs with values of N that violate Eq. 4.28. Note that in both the (5 = 1,N = 1203) T A S E P and the (5 = 2,N = 1204) T A S E P the last defect is located at lattice site — 2 and the TASEPs share a common value for 13*. Additionally, TASEPs (6 = %N = 1203), (6 = 3>N = 1204), and (8 — 1, N — 1205) all share a common location for the last defect site, and share a common value of j3*. 2 Chapter 4. A Totally Asymmetric Exclusion Process with Periodic Structure 61 Fig. 4.12 show that violations of condition 4.28 produce a phase diagram that does not preserve the symmetry about the a = 8 line that is found in the standard T A S E P phase diagram. Conversely, when Eq. 4.28 is satisfied, this symmetry re-appears. Note that, within simulation accuracy, the phase diagrams of the non-symmetric dual-rate TASEPs retain the general three phase form the standard T A S E P phase diagram for all values of 8 and N tested. The values of a* and 8* for the non-symmetric phase diagrams of Fig. 4.12 suggest that the location of the first defect site determines the value of a*, while the location of the last defect site independently determines the value of 8*. For example, consider the T = 3 case with (8 = 1, N = 1203) where a* « 0.27 and B* « 0.19. By construction, this realization of the T = 3 dual-rate T A S E P and the symmetric version of the T = 3 dual-rate T A S E P with (5 = 1,N = 1205) both place the first defect at lattice site 1. The two realizations of the T A S E P then share the same value for a* (~ 0.27). Conversely, the non-symmetric (5 = 1, N = 1203) case and the symmetric (5 = 2, N = 1204) case both have the last defect in lattice site N — 2, and share the same value of 8* 0.19. This pattern held for all values of 8 and N examined. Finally note that the value of a* for the 8—1 T A S E P is greater than the value of p 2 [a* « 0.27, p = 0.25). When we consider 2 that the first lattice site of the T A S E P is a defect site, then the effective injection rate into the remainder of the T A S E P lattice cannot exceed the defect rate. As a result the T A S E P will be in an entry limited phase until a at least matches p . 2 4.6 Summary We have developed approximate, mean field, methods to compute the steady-state current and densities of totally asymmetric exclusion processes involving two internal hopping rates. Additionally, we have simulated the two-rate T A S E P and have explored its phase diagram. We find that the dual-rate T A S E P retains the three phases found in the standard T A S E P model. Within each of these phases the best of our mean field theories provide reasonable approximations to the particle currents and densities. In particular, a maximal current phase mean field theory was developed that provides accurate estimates Chapter 4. A Totally Asymmetric Exclusion Process with Periodic Structure 62 (Eq 4.14) for the current and the density (Figs 4.4(a) and 4.7) in the (5 = 2 , T = 2) case. Similarly, in the boundary limited phases a mean field method was introduced for the (5 = 2 , T = 2) case (Eqs 4.20 and 4.21), which produced accurate estimates for the boundary currents. Additionally, we have developed a primarily numerical mean field approach that we have termed the finite segment mean field method (FSMF). This method can rapidly produce estimates for the currents and densities in all three phases of the dual-rate TASEP. However for larger periods, the F S M F method was found to produce increasingly inaccurate current estimates. This suggests the method may be most useful for treating TASEPs with highly localized regions of non-uniform hopping rates, rather than the periodic case investigated here. Finally, we have used Monte Carlo simulations to briefly investigate TASEPs with alternative arrangements of periodic defects. Specifically, we investigated dual-rate TASEPs where the first defect was located at an arbitrary position 5 within the first T sites of the lattice. We found that the dual-rate T A S E P possesses particle-hole symmetry when the number of sites in the T A S E P lattice equals mT + 25 for an integer m > 0. When the size condition was not met, the phase diagram of the dual rate T A S E P was no longer symmetric about the a = B line. Additionally, results from Monte Carlo simulations suggest that in the non-symmetric case, the value of a* is determined by the position of the first defect in the T A S E P lattice, while the value of B* is independently determined by the position of the last defect. 63 Chapter 5 Clustered Bottlenecks in m R N A Translation and Protein Synthesis 5.1 1 Background Having developed techniques for analyzing TASEPs with spatially varying hopping rates in Chapter 4, we now apply a variation of the dual rate T A S E P to the topic of biological protein synthesis. During protein synthesis, ribosome molecules bind at one end (the 5' end) of a messenger R N A (mRNA), translate along the m R N A sequence, and detach from the mRNA with the completed protein product at the 3' termination end (Fig. 5.1). Each translation step requires reading (translating) a nucleotide triplet (a codon) and the binding of a freely diffusing transfer R N A (tRNA) molecule carrying the amino acid complementary to that codon [7]. Besides being a critical final stage of gene expression, control of protein synthesis is vital for protein adaptation and evolution [73-75], for the control of viral parasitism [76], and for high yield, cell-free, synthetic protein production [77]. As the correct transfer R N A must bind to the ribosome m R N A complex for a codon to be successfully translated, the overall rate of translation is strongly dependent on the concentration of transfer RNAs in the cellular cytoplasm. For example, so-called "slow" or "defect"codons (those with rare corresponding tRNA and/or amino acids), are known to inhibit protein production [78, 79]. These defect codons typically include C T A (Leucine), A T A (Isoleucine), A C A (Threonine), C C T and C C C (Proline), C G G , A G A , and A G G (Arginine) and arise infrequently (about 3-4% of mRNA codons in E. coll) [80, 81]. The concentration of the t R N A molecules complementary to these defect codons can be as much as 30 times lower than the concentration of the more abundant tRNA species [8]. T h e work in this chapter is based on "Clustered Bottlenecks in m R N A Translation and Protein Synthesis",Tom Chou, Greg Lakatos, Phys. Rev. Lett. 93 (2004) 198101 x Chapter 5. C l u s t e r e d B o t t l e n e c k s i n m R N A T r a n s l a t i o n a n d P r o t e i n Synthesis 64 Figure 5.1: The mRNA translation/protein synthesis process. Ribosomes move unidirectionally along the mRNA as tRNAs (not shown) deliver the appropriate amino acid to the growing protein. Codons with low concentrations of corresponding t R N A ("slow-codons") result in bottlenecks that locally slow ribosome motion. A totally asymmetric simple exclusion process on N lattice sites is used to model the mRNA translation. Slow-codons are represented by lattice sites with a reduced particle movement rate p < pi, 2 Chapter 5. Clustered Bottlenecks in mRNA Translation and Protein Synthesis 65 Defect codons can appear throughout the mRNA, be at or near the initiation site, the termination site, or in the interior region [69]. of rare codons near the 5' binding site of E. coli Statistics indicate a higher occurrence genes [82]. Even more striking is the tendency for rare codons to cluster [80]. Clusters of approximately two to five consecutive defect codons occur frequently in E. coli, D r o s o p h i l a , yeast, and primates [80]. The existence of defect codons suggests at least two possible mechanisms for regulating protein production. The first regulatory mechanism involves the dynamic increase or decrease in the concentration of specific tRNAs in a cell, in response to a changing environment. As the concentration of a particular tRNA dropped, the rate of translation of any mRNA rich in the codon complementary to the targeted t R N A would also decrease. The second regulatory pathway is more evolutionary in nature, and would be useful to an organism like a virus that has no transcriptional machinery of its own. In the absence of such machinery, a virus can regulate the production of its proteins by adapting its genetic code to use codons complementary to tRNAs that are more or less abundant in the host organism [83]. This type of adaptation would be expected in transgenic viruses that jump between different host species. Finally, we note that an understanding of the effects of varying tRNA concentrations has also been useful in applications of recombinant D N A technology. For example, the genetic code for green fluorescent protein (naturally found in various species of jellyfish), has been modified to use codons with conjugate tRNAs that are in high concentration in mammalian cells. This modification increased the utility of recombinant green fluorescent protein in several experimental studies of mammalian cell lines [84]. 5.2 Model Although the strength, number, and positioning of defects can affect local ribosome densities and overall translation rates, to date there has been no quantitative model describing how various defect sequences can control ribosome throughput and protein synthesis. Here, we use the TASEP, continuous time Monte Carlo simulations based on the B K L algorithm [67], and the finite-segment mean field method to model how specific codon us- Chapter 5. C l u s t e r e d B o t t l e n e c k s in m R N A T r a n s l a t i o n a n d P r o t e i n Synthesis (a) Pi 7> (b) • L = 5 k= 1 66 ^ L = 77. £ = 5 Figure 5.2: Placements of slow-codons or "defects" (blue shaded lattice sites), (a) Two adjacent defects subsumed in a L = 5 finite segment, (b) Two defects separated by four normal codons and subsumed in L = 11 finite segment. ages that give rise to local delays in translation can be used to suppress protein synthesis. In mapping the T A S E P onto the mRNA translation process, the T A S E P lattice represents the mRNA, and each T A S E P particle represents a translocating ribosome. Each time a T A S E P particle makes a hop, a codon is read, and a tRNA delivers its amino acid to the growing polypeptide. When a T A S E P particle reaches the right side of the lattice it detaches at a rate 8 as normal, and a completed protein would theoretically be available for use in the cell. Note that mRNAs are typically 100 to 1000 codons in length, leading to T A S E P lattices approximately the same size as those investigated in the preceding two chapters. The specific variant of the T A S E P used to model a defect bearing mRNA is displayed in Fig. 5.1. As defect codons are relatively rare in typical mRNAs (with a probability of occurrence of approximately 0.03/codon) and typically cluster, the model shown in Fig. 5.1 only exhibits a finite number of localized defects, rather than the periodic arrangement considered in Chapter 4. The nomenclature used to describe the spacing of defect codons and the size of the finite segment used in the mean field computations that follow is displayed in Fig. 5.2. As we work with finite, localized clusters of defects rather than periodically arranged defects, the application of the F S M F method to the ribosome problem differs somewhat from the methodology used in Chapter 4. To see how the F S M F method is used in this application consider the interior segment of the T A S E P lattice displayed in Fig. 5.2(b). We define the boundaries of the defect cluster to be marked by the defect that is closest to the a injection site on the left, and the defect that is closest to the B extraction site on the right. If we now assume that far from the defect cluster the T A S E P will behave like a standard defect free TASEP, we can defineCT_to be the density in the region of uniform Chapter 5. Clustered Bottlenecks in mRNA Translation and Protein Synthesis hopping rates to the left of the defect cluster. Similarly we define o + 67 to be the particle density in the region of uniform hopping rates to the right of the defect cluster. If we then apply the simple mean field approach that was successful in the standard T A S E P to the regions to the left and right of the defect region, we find J = cx_(l — cr_) = o (l + — o ), + where the hopping rate i in the uniform regions has been normalized to one. This implies P or cr_ = (1 — o~ ). Given that the T A S E P lattice contains defects, we that either cr_ = o + would expect that o + + < cr_ and we rule out the first of these two solutions as unphysical. If we now encompass the defect cluster with a finite-segment that extends into the uniform hopping rate regions found to either side of the cluster, we find a ff = (3 s = (1 — o + ). e e Recall that in the standard TASEP, far from the boundaries the occupancy correlations between lattice sites decay to zero. If we assume that this holds for the regions far from the defect cluster then the mean field coupling of the finite segment to the rest of the T A S E P lattice should be an extremely good approximation. With these assumptions, the F S M F method proceeds as described in Chapter 4, but here we only need to find a fixed point in the single density a . + Thus far we have only described the application of the F S M F method to the case where the defect codons in the mRNA are far from the mRNA ends. However, if we assume that far from the defect cluster the T A S E P behaves like a standard defect-free T A S E P , then lattices with localized defects near the injection or extraction ends would behave qualitatively like standard TASEPs with reduced values of a or (3 respectively. As a result, we focus on defects located far from the boundaries. Additionally, in a typical E. coli bacterium, there are roughly three times as many ribosomes as there are t R N A molecules of the most abundant variety [8]. This suggests that the movement of ribosomes rather than ribosome binding and release, is the rate limiting step in the translation process. Consequently, we focus our analysis on TASEPs with large values of a and ft, that in the absence of defects, would be in the maximal current phase. Chapter 5. C l u s t e r e d B o t t l e n e c k s in m R N A T r a n s l a t i o n and P r o t e i n Synthesis i 0.3 0.25 ^ 0.2 ^ 0.15 i i i | i i | i 1 1 1 | 1 I I 1 | 1 T [ [ (a) Current, single defect i i 0.1 0.05 0 1 / , , , , ;\ I , , , , 0.2 0 ,1 6 8 I , , , , 0.4 I , , • - • M C Simulation - - L = 1 FSMF L = 3 FSMF - L = 5 FSMF 0.6 0.8 1 1 1 1 1 1 1 1 1 1 1 1 (b) Current, m contiguous defects i 1 = 0.01 — P == 0 . 1 0 • - * p == 0 . 2 0 2 cTO.75 3* 2 0 .1•—.5 , - _ _ _ _ J ^ : ~ P _-AP ^^^^ ~ *• 0.25 ' 1 " 1 t 2 2 2 == 0 . 3 0 == 0 . 5 0 l t 1 l 1 1 l l l l 1 3 4 l i i 5 m Figure 5.3: (a) Comparison of steady-state currents Ji(p ) for a single defect derived from M C simulations with those obtained from a //-segment F S M F method. For L > 5, the F S M F results are within 2% of those from M C simulations. The boundary injection and extraction rates were not rate limiting (a = 6 = 10). (b) Further reduction of the steady-state current as successive, identical defects are added. The first few defects cause most of the reduction in current. Note that the ribosome current varies significantly as ~> is changed. Thus for clarity, the currents in the figure have been normalized by the slow-codon hopping rate (p ). Also note that for small ~>, J / p ^ 1/2 consistent with Eq. 5.3. 2 2 2 2 2 2 Chapter 5.3 5. Clustered Bottlenecks in mRNA Translation and Protein Synthesis 69 Results and Discussion We begin by considering the effect of a cluster of close-packed defect codons in the centre of an mRNA on the overall rate of protein production. Intuitively we would expect that the current J will diminish upon addition of successive, contiguous defects. To determine the approximate functional form of the current reduction, we can work out the finitesegment mean field approximation to the ribosome current to lowest order in the defect codon translation rate p . 2 Let m be the number of contiguous defects and as before let cr equal the density downstream from the last defect. If we set the finite segment to + exactly cover the m contiguous defect codons we can set the effective injection rate to eff = (Pi/p )(l a 2 — cr+) and the effective extraction rate to /? ff = (1 — c ) . e + Making use of the exact matrix product solution from [15] for a T A S E P of length m with an internal hopping rate of 1, the finite segment expression for J is J Where R m j, ' _ R -i(P~j) - Rm-i(u~i) ~ D (ro-l\ c / -IN r m!(m r m+l . 1 = E — M l v lJ-, + 1 _ } l/a ff e -„ n + 1) (n — .l)(2m. — .n)\,- n=2 ( . \P ) A m+l t— + • is given in Eq. 2.2. Now consider the terms in J that depend on = As o , m P 2 l ) . ! f .„. e / a l m " - > + (5.2) , v I ° ' ° ' » ' ' is the density downstream from the block of contiguous defects, in the limit of small p , o+ must approach 0. Thus the lowest order term in the power series expansion 2 of <7+ must be at least linear in p . The sum in Eq. 5.2 starts from 2 and thus 2 R (l/a g) m e and Z? _i(l/a fj) are both of at least order 2 in p and can be ignored in the computation m e 2 of the lowest order term in J. Then, lowest order in p 2 Chapter 5. Clustered Bottlenecks in mRNA Translation and Protein Synthesis 70 ( n - l ) ( 2 m - (n + 2)! i^-!^ ) 1 E ^ ^U&ff) ( -l)!(m-n m ^ n=2 (n-l)(2m-rc)! fm)!(m + 1 — ri)\ (5.3) m + l As expected in the limit p —* 0, J approaches 0. Similarly, in the limit of large m, the 2 current predicted by Eq. 5.3 approaches J = P 2 / 4 , the current of a long uniform lattice with hopping rates p in the maximal current regime. Figure 5.3(b) plots J (p )/p 2 m 2 2 for various p as a function of m defects using a large L = 14 finite segment. Note that as 2 the currents for a wide range of p values are displayed, all currents in Fig. 5.3(b) are 2 normalized by the slow-codon hopping rate p . For all p < 1, J (P2)/P2 2 2 m approaches 1/4 as m —> N. For strong defect codons (p < 0.3), the largest decrease in J {p ) 2 m occurs as m goes 2 from one to two. Thus, we consider the effects of placing only two defects in the mRNA interior. Figure 5.4(a) shows the expected currents J (p ,k) 2 2 across a lattice containing two defects spaced k sites apart. Results for k < 5 were computed using an L — 14 F S M F method, while those for k > 5 were obtained from M C simulations. The largest reduction in the current occurs when two defects are spaced as closely as possible. Considering the ribosome densities near the defects, we find that as k increases, the density downstream from the first defect recovers to the bulk value cr_ before encountering the second defect. This behaviour is clearly shown using both F S M F and M C simulations in Fig. 5.4(b) for a pair of defects (p = 0.15). Also note that Fig. 5.4(b) shows that the density profile 2 between the two defects is approximately linear, similar to the known behaviour of the standard T A S E P along the a = (3 line in the phase plane. As before, we can analyze the finite segment mean field approximation to the current in the case of two defect codons separated by k sites. Again define o + to be the density to the right of the second defect, and consider a finite segment of k sites that covers the region that immediately follows the leftmost of the two defective sites, and extends up to and includes the second defect site. Normalizing the rate of transcription of the normal Chapter 5. Clustered Bottlenecks in mRNA Translation and Protein Synthesis 71 codons to one, the effective injection and extraction rates are a ff = ArfF = P ( l ~ "+)• e 2 Applying the exact matrix-product method solution from [15] to the finite segment gives -1M J = lim R (8-^) - k (5.4) R (a-J) k As the limit in Eq. 5.4 is an indeterminate form, we apply l'Hopital's rule to find ™ &=l/<*ef In Eq. 5.5 7 (/c) = ^~fc+i ^)^ • As o;fr ~ n 2 ! n e P . 2 /c+i (5.5) (fc+l)-n y^7n(fc)^ ff n=2 e 2 the lowest order term in the numerator of Eq. 5.5 must by linear in p , while the lowest order term in the denominator must be 2 independent of p . So that, to lowest order in p , 2 2 j ^ j (k k - l)o! lk+i{k) eff _ fc kp (l 2 - o) + + l (5.6) _ kp 2 fc+l' In going from the first to second line in Eq. 5.6 we could drop the o + term in the numerator as o , the density downstream of the second defect, must be at least order one + in p . 2 Equation 5.6 shows that the current for two identical, strong, defects is at most a factor of 2 smaller than that for a single defect. This maximal difference occurs when P 2 —» 0 and when the two defects are adjacent to each other. variation in J2GD2, k)/J\(p ) 2 Figure 5.5(a) plots the as a function of k for various p , and shows that in the limit 2 of large k the current with two separate defects approaches the current produced in the presence of a single defect. Also note that the non-trivial defect spacing dependence of J is a consequence of the hard-core exclusion between T A S E P particles. If the particles had no exclusion and were free to pass by each other, their complete traversal time would be the sum of the individual hopping times at each site, and would not depend on the defect arrangement. Chapter 5. C l u s t e r e d B o t t l e n e c k s i n m R N A T r a n s l a t i o n a n d P r o t e i n Synthesis 72 i Figure 5.4: (a) Steady-state currents across a lattice with two identical interior defects spaced k sites apart. The current is suppressed most when the defects are closely spaced, (b) The density profiles near a pair of defects (p = 0.15) of various spacings k from F S M F computations. The circled points correspond to the defect positions for k = 6. For larger k, density boundary layers decay, allowing particle build-up behind the second barrier, enhancing the current. 2 Chapter 5. Clustered Bottlenecks in mRNA 1 0 5 10 1 15 Translation 1 20 m 1 25 and Protein Synthesis 1 73 1 30 35 40 Figure 5.5: (a) The dependence of the normalized steady-state current Jiipi, k)/J\(p2) as a function of the defect separation, (b). The upper bound for the mean current of a lattice with m defects randomly distributed within the central N = 3 0 0 sites. A lower bound as m —> N is J(m —> A Q R A N — ^ 2 / 4 . Chapter 5. Clustered Bottlenecks in mRNA Translation and Protein Synthesis 74 From the two-defect current J , we can estimate the expected current for m randomly 2 distributed defects. Using Eq. 3.11 from Chapter 3, the total number of ways m defects can be placed on N interior sites, such that the minimum pair spacing is greater than or equal to k is ~ / (N -(m \ »rs Z (m,N)=l The 1)\ ). J JK 5.7 J probability that at least one pair of defects is separated by exactly k sites is then Q (m,N)= { k As - l)(km 1 k Z k ^i Z k + l ) . (5.8) the ribosome current in an mRNA with two or more defect codons all separated by at least k normal codons must be less than or equal to J2{P2, k), we find the following upper bound on the current L(N-l)/(m-l)J J(m,N,p ) K 2 RA < Qk{m,N)J (p ,k). 2 (5.9) 2 k=l This upper bound will be more accurate if the defect density is low enough that the probability of the defects forming closely-packed clusters is very small. J(m,N = 300)RAN normalized by J\{p ) 2 is shown in Fig. 5.5(b). the current is approximately that of a single defect. The current For very small m, The current is most sensitive to the number of random defects when m is approximately 10, which corresponds to m/N approximately equal to 0.03, roughly the fraction of slow codons observed in mRNAs of E. coli. However, the observation of nonrandom clustering in real mRNAs suggests that using the typical biological distribution of defect codons would yield ribosome currents well below the upper bound given by Eq. 5.9. 5.4 We Summary have found that not only can a single defect directly inhibit translation across it, but a few defects properly distributed, can further slow protein production by a factor of approximately two to four. Although maximal current reduction is achieved by clustering Chapter 5. Clustered Bottlenecks in mRNA Translation and Protein Synthesis 75 defects as closely as possible, successive additions beyond approximately five defects does little to further reduce the current. Additionally, defect codons which are all spaced more than approximately eight sites apart will not reduce the ribosome current significantly more than a single defect codon. Our results are qualitatively consistent with the idea that a single, localized region with defect codons, becomes the rate limiting step for mRNA translation. The existence of slow codons near the mRNA start codon [82] suggests that their role is to suppress protein synthesis. Conversely, we speculate that an mRNA with a finite number of well-separated defects located at appropriate positions in the interior of the mRNA, may provide pause points for local successive protein folding with minimal reduction in the overall ribosome current. 76 Chapter 6 First Passage Times of Driven D N A Hairpin Unzipping 1 6.1 Background Over the last two decades, D N A biotechnology and nanofabrication have been transformed by advances in single molecule manipulation techniques. In particular, there has been strong interest in the mechanisms of single stranded D N A transport across membranes [37-41], as well as investigations of force-pulling mechanisms for inducing D N A dehybridization ([48-50]). Here, dehybridization is defined as the separation of the two stands of double-stranded D N A , accomplished by breaking the bonds between complimentary D N A bases. While the process of single stranded D N A transport through membrane channels has been extensively studied [44-46, 49], only recently has the process of driven double stranded D N A dehybridization by transmembrane nanopores been investigated [51-53, 85]. In this chapter, the dehybridization of D N A resulting from the pulling of D N A molecules through a transmembrane channel is studied. This work is primarily motivated by the measurements of Mathe et al. [52], and related to the experiments of Nakane et al. [53], and Sauer-Budge et al. [85]. In the experiments of Mathe et al, a transmembrane potential was applied to drive the translocation of D N A hairpins. As the transmembrane voltage was applied, the ionic current through a trans-membrane a-hemolysin channel was simultaneously monitored. When any part of the D N A strand was inside the transmembrane channel, the channel was blocked, and the ionic current fell below the open channel value. When the D N A strand exited the channel, open-channel ionic currents resumed and were measured. The experiments measured the distribution of first passage The work in this chapter is based on "First Passage times of driven D N A hairpin unzipping", Greg Lakatos, T o m Chou, Birger Bergersen, and Gren N . Patey, Phys. Biol. 2 (2005) 1-9 l Chapter 6. First Passage Times of Driven DNA Hairpin Unzipping 77 times for the escape of the D N A from the a-hemolysin channel. The mean first-passage times r(V) were then derived as a function of driving voltage V. The driving voltage- dependent mean first-passage times of perfectly complementary hairpins were compared to those of D N A hairpins with single defects (non-homologous base pairs). Since there are effectively fewer base pairs to break, D N A hairpins with non-complementary defects had a shorter first-passage time. The first-passage time data from these experiments [52, 53] have all been analyzed by modeling the D N A dehybridization as an activated process with a single free energy barrier. In these models there are three parameters: an attempt frequency (prefactor), a barrier free energy, and a partial charge per D N A base. The energy barrier should be approximately the difference in free energy between hybridized and fully dehybridized DNA. The D N A hybridization energetics can be modeled using both a simple free energy model where each type of hybridized base (i.e. guanine-cytosine (GC) and adenine-thymine (AT) base pairs) contributes a constant amount to the total free energy, and using free energies produced by more sophisticated hybridization models such as the Multiple-Fold (MFOLD) 2-state hybridization server [4]. Permitting a 1.5 to 2 order of magnitude variation in the attempt frequencies [52], the single barrier fits give a good match to the data. In order to understand the first passage times observed in these experiments and investigate possible deviations from log-linear behaviour at high and low voltages, we obtained estimates for the mean voltage-dependent first passage time, r(V), using a twodimensional multi-barrier stochastic model. The D N A hybridization energetics are modeled using both a simple free energy model where each G C base pair contributes 3k T, and B each A T base pair contributes 2/CBT ( T = 288K) to the total hybridization free energy, as well as with the free energy model used by the M F O L D 2-state hybridization server [4]. The predicted mean first passage times show the same nearly linear dependence of log(-r(V)) on the applied voltage at moderate voltages ( « 25 to 250mV), as was observed by Mathe et al. for voltages between 30 and 150mV. However the model predicts an additional behaviour in the low voltage limit. For voltages below approximately 25mV, our model suggests that in the experimental setup of Mathe et al, sliding of D N A in the reverse direction (the direction opposite to the applied electric field) is the dominant Chapter 6. First Passage Times of Driven DNA Hairpin Unzipping 78 mechanism of D N A escape. 6.2 Model We model the process of D N A extraction through a membrane pore using the zipper model displayed in Figs. 6.1 and 6.2. The model assumes that the D N A forms a hairpin consisting of base pairs in the self-hybridizing region, and that the hairpin is connected to a single-stranded poly-A tail NA bases in length (see Table 6.1 for the definition of NA and the other size parameters used in the model). In the experiments of Mathe et al. N = 10, and NA = 50. The dehybridization and translocation states of the model are labelled with two indices, no and n\. If the bases of the D N A strand are enumerated starting from the 3' end, index n indicates the position of the D N A strand relative to the 0 start of the transmembrane region of the cv-hemolysin pore. As the D N A strand can slide in reverse, no can take on negative values. Similarly, index n equals the number of hairpin x base pairs that have been dehybridized (See Fig. 6.1(b)). Using this nomenclature, the D N A hairpin is fully dehybridized when n\ = N. As the dehybridized D N A hairpin can escape the pore region by sliding in reverse (in the cis direction, see Fig. 6.1(b)), we make all states with n = —NA absorbing. Forward escape (escape into the trans chamber) can 0 only occur after the D N A hairpin has completely dehybridized, and so we also make state (n = 2N 0 + N LOOP + N, L Defining P(w,v,t) n — x N) absorbing. to be the probability of being in state (n = w,ni 0 = v) at time t, the system then evolves according to a master equation dP(w, v, t) ~ ^L(i,j\w,v)P(i,j,t) dt In Eq. 6.1 the expression L(i,j\w,v) - P(w,v,t)^2L(w,v\i,j). (6.1) is the rate of transition from the state with (n = 0 i^m = j) to the state (n = w,ni = v). We restrict the set of allowable state transitions 0 to include only four moves: sliding of the entire D N A strand in the trans direction (no —> (n + 1)), sliding in the cis direction (n —> (n - 1)), dehybridizing a D N A base pair 0 0 0 (ni —» (ni +1)), or rehybridizing a D N A base pair (m — > {n - 1)). In all cases we assume x Chapter 6. F i r s t Passage T i m e s of D r i v e n D N A H a i r p i n Unzipping Parameter Description V Number of base pairs in the D N A hairpin Number of bases in the poly-A tail attached to the 3' end of the D N A hairpin (50). Number of bases that can simultaneously occupy the transmembrane part of the pore (12) Number of bases that can simultaneously occupy the vestibule part of the pore (16) Number of bases in the hairpin loop (6) N A N L N v A^loop MR The 79 rate of rehybridization of a completely denatured hairpin (set to 0) Table 6.1: The fixed size and rate parameters used in the two dimensional stochastic model. The size parameters are set by the structure of the a-hemolysin pore and the D N A strand, and are not used as fitting parameters. Figure 6.1: (a) Division of the a-hemolysin channel into the cis, vestibule, transmembrane, and trans regions. The D N A hairpin is inserted into the channel from the cis region. The number of D N A bases that can fit into the vestibule region is set to N = 16, while the number of D N A bases than can fit into the transmembrane region is NL = 12. (b) A n example of the D N A hairpin translocation experiment. The process of D N A translocation is described by two coordinates no and n\. no gives the position of the first hairpin base relative to the mouth of the transmembrane part of the channel, while index n-[ gives the number of hairpin base pairs that have been dehybridized. V Chapter 6. First Passage Times of Driven DNA Hairpin Unzipping 80 CD C A A A 7 6 5 A 4 A A 3 A 2 A A A A 3 ' 1 Base Index Figure 6.2: The enumeration scheme for D N A base pairs in the hairpin and the associated state space for the 2D stochastic model. In this example N = 7, N = 4, and there is a defect at base 4. States with rt\ = 3 are inaccessible as the presence of a defect at site 4 makes it impossible to have exactly 3 base pairs dehybridized. Transitions between states with different values of ri\ have attempt frequency u. In all cases, transitions between states with different values of n have attempt frequency p.. States with n — 2N + N\ + N (a D N A strand pulled completely into the trans chamber), and n = —NA (a D N A strand pulled completely into the cis chamber) are absorbing. A 0 Q OOP 0 L Chapter that L(i,j\w,v) 6. First Passage Times of Driven DNA Hairpin Unzipping 81 has the form L(i,j\w,v) = f(iJ\w,v)e-^ m m ^ ^ G w v ) m (6.2) . Equation 6.2 has the standard form of a rate in an activated process [86] where f(i,j\w,v) represents the attempt frequency for escape from state while AG(i, j\w,v) fine AG(i, j\w,v) to state is the free energy difference between states (w,v) and (w,v), To de- we divide the pore-membrane complex into four regions: the cis region, the vestibule region, the transmembrane region, and the trans region (see Fig. 6.1(a)). Furthermore, we assume that the total free energy of the D N A strand is the sum of the contributions from the parts of the strand in each of these four regions. We begin by considering the major entropic contributions to the total free energy of the D N A strand, from each of the four regions. In the trans region, all the D N A is single stranded and we set its contribution to the free energy equal to that of a polymer tethered to the membrane ([42, 44, 47]). The contribution from the D N A bases in the trans region to the change in free energy upon going from state (i,j) to state (w,v) is then (6.3) In Eq. 6.3, N (no,ni) trans is the number of D N A bases in the trans chamber in state (n ,ni), and b is the persistence length of single stranded D N A (approximately two inter0 base distances) [49]. Since the bases in the transmembrane region of the pore are highly confined, there will be an entropic penalty for storing bases in this region. We assume this entropic penalty is ~ lkB per base [47]. The contribution to the total free energy change from the bases in the transmembrane part (TMP) of the pore is then A G M P ( « , j\w,v) T where A " 7 TMP = k T(N (w,v) B TMP - N (i,j)) TMP (6.4) ( n , n\) is the number of single-stranded D N A (ssDNA) bases stored in the 0 Chapter 6. First Passage Times of Driven DNA Hairpin transmembrane part of the pore in state ( n , n i ) . Unzipping 82 Turning to the vestibule region we 0 note that the diameter of the vestibule varies from approximately 26A at the vestibule entrance, to 36A near the middle of the vestibule region [87]. Though the vestibule region is less confining than the transmembrane part of the pore, we nonetheless expect that storing ssDNA bases in the vestibule will carry an entropic penalty. The contribution to the total free energy change from ssDNA bases in the vestibule region is AG (iJ\w,v) = s T(N (w,v) vest where A\/ tibuie(^o, \ ) n es v - N (iJ)), vest (6.5) vest is the number of ssDNA bases located in the vestibule region, including all the single stranded bases created when part of the hairpin dehybridizes. The parameter s v > 0 is the entropic cost for storing a D N A base in the vestibule region. While the hairpin is in the vestibule region and s v is non-zero, the free energy cost to open a base pair is increased by 2s T and the hairpin is stabilized. While we leave s as a v v fitting parameter, the confinement free energy can be estimated from the entropy cost per persistence length to confine a polymer to a narrow tube: s ~ k,B{l/D)i v [88]. Using the mean vestibule diameter of D = 3.1nm and a persistence length / = 0.8nm, the entropic cost of confining a base to the vestibule is approximately 0.1/cg. Next, we consider the cis region, where we treat the single stranded D N A as a polymer tethered to the membrane. With this free energy, the contribution to A G from the cis region is ([42, 44, 47]) AG (i,j\w,v) cis where N (n , cis 0 = ^k Tln B N (w,v) cis (6.6) Nds(i,j) rii) is the number of single stranded bases located in the cis region, exclud- ing the single stranded bases that are produced by partial dehybridization of the D N A hairpin. Next, we introduce the enthalpic contribution to the free energy from the change in the electrostatic potential associated with sliding a charged D N A strand AG (i,j\w, el v) = U(w,v, V) - U(i,j, V). (6.7) Chapter Here U(n , n 0 X) 6. First Passage Times of Driven DNA Hairpin Unzipping 83 V) is the potential energy of the D N A strand and is given by N +n A 0 £ U(n , ,V)= 0 ni (6.1 qcP(i,V) i=l b where (fr(i, V) is defined to be the potential at a distance i inter-base spacings from the cis entrance to the transmembrane region, q is the mean partial charge per D N A base, and l 0 is the location of the single stranded base that is closest to the cis chamber. We assume that the entire potential drop V occurs across the w i d t h of the transmembrane part of the channel [44], and that the voltage profile i n the transmembrane channel is linear 0 1/2)^V \ if i < 0 ifO<i<N (6.9) L if i > N L To complete A G we include the free energy cost associated w i t h dehybridizing a D N A base. We determine this cost using two different models: the first is a simple free energy model where the binding free energy of each A T and G G base pair is set to 2kBT 3kBT, and respectively (the "2-3" free energy). Alternatively, we also use the free energies predicted by the M F O L D D N A folding program [4]. W i t h either of these models and Eqs. 6.3-6.7, the total change i n the free energy resulting from a transition from state (i,j) to state (w, v) is AG(i,j\w,v) = (g(v) - g(j)) + AG (i,j\w,v) cis AG p(i,j\w,v) TM AG (i,j\w,v)+ e] + AG (i,j\w,v) vest + AG trans + (6.10) (hj\w,v). Here g(ni) is the free energy of a hairpin in bulk solution w i t h n\ bases dehybridized, and is given by either the M F O L D program or the "2-3" free energy model described above. Having defined the energetics of the state transitions, we are left w i t h the attempt Chapter 6. First Passage Times of Driven DNA Hairpin Unzipping 84 frequency in Eq. 6.2. The model has two basic attempt frequencies associated with the two types of transitions: p, the attempt frequency for sliding transitions (no —» no ± 1), and u, the attempt frequency for hybridization transitions (n\ —> ri\ ± 1) (see Fig. 6.2). To develop a rough estimate for p, we set this attempt frequency equal to the inverse of the relaxation time of the single stranded D N A confined to the transmembrane part of the channel. The attempt frequency is then [47] ^CJnTC ( 6 - N ) where a is the inter-base spacing of single stranded D N A , 7 is a measure of the stiffness of the D N A in the channel, and Q is the typical friction coefficient for a base in the channel. Using the mean values for £ and 7 given in [47], we find p ~ 10 s . Using the 5 _1 value of C given in [44] for a persistence length of single stranded D N A confined in the transmembrane region gives an estimate of p ~ 10 s~" . We constrain p to lie within the 6 x range of these two estimates. To find an estimate for u>, we use a Morse-like potential to model the hydrogen bonds between two D N A bases [48]. A calculation of the Kramer's escape rate to open a single base then gives a rough estimate of u> » 1 x 1 0 s , consistent 12 _1 with previous studies of D N A dehybridization [48]. While the D N A extraction process is not an equilibrium process, and violations of detailed balance can affect the mean first passage time appreciably, we assume that the processes is near enough to equilibrium for detailed balance to apply. The one exception is the rate of rehybridization of a completely dehybridized hairpin. In this case we assume that once the D N A strand has completely dehybridized it escapes the pore complex before rehybridization can occur, and set the rehybridization rate of the last hairpin base to 0. Having completely defined a two dimensional stochastic model for the hairpin extraction process, we are left with four free parameters summarized in Table 6.2, which we fit to experiments. This model does not include a number of features of the double stranded D N A translocation and dehybridization process. For example, by assuming the D N A hairpin dehybridizes only from the free end, the effects of bubble formation are ignored by our model. Here, bubbles refer to dehybridized regions of D N A that are bounded by at least one Chapter 6. First Passage Times of Driven DNA Hairpin Parameter Unzipping Description A t t e m p t frequency for sliding transitions (n —> n ± 1) A t t e m p t frequency for hybridization transitions (rii —> n i ± 1) Mean partial charge per D N A base Entropic penalty for confining a s s D N A base to the vestibule region 0 u q Sy 85 0 Table 6.2: List of fitting parameters used i n the two dimensional stochastic model. hybridized base on b o t h sides. Under physiologic conditions, structural fluctuations i n double stranded D N A lead to the formation of bubbles that are typically on the order of several tens of base pairs [89]. Our treatment and the measurements of Mathe et al consider the dehybridization and extraction of D N A hairpins w i t h base-pairing regions approximately 10 bases in length, and thus considerably smaller than the typical bubble size found i n D N A at experimental conditions. Additionally, recent experiments [90] have studied the formation of bubbles in small, atypically bubble-prone p o l y - A T D N A molecules approximately 18 bases in length. Here, the formation of D N A bubbles approximately 2 — 10 bases i n length was observed i n the D N A strand. However, Altan-Bonnet et al. [90] were able to measure the typical lifetime of these small bubbles to be approximately 50^,s. T h i s lifetime is typically much smaller than the mean time for the dehybridization and extraction of small D N A strands through a-hemolysin at moderate driving voltages. Thus we assume the effects of bubbles can be subsumed into the dehybridization attempt frequency to. Additionally, the model does not include partial-registry binding, where the hybridization region of the hairpin is offset by one or more bases (i.e., i n the arrangement of F i g . 6.2, complement base 2 hybridizes w i t h base 5). Under the assumption that the doublestranded D N A complex is well equilibrated prior to the start of the extraction processes such out of registry binding is unlikely, due to its high energetic cost. Furthermore, in this model of D N A dehybridization and sliding, we have assumed that the dominant enthalpic contributions to the free energy of the D N A strand come from the D N A hybridization, and the applied voltage. In doing so, we have ignored the possible interactions between the D N A stand and the lipid membrane or parts of the a-hemolysin channel. Finally, previous studies of D N A translocation [44, 91] have shown that the use of an effective free energy of the kind described by E q . 6.10 is strictly valid only when the Chapter 6. First Passage Times of Driven DNA Hairpin Unzipping 86 relaxation time of the D N A strand is much shorter than the typical time required for the strand to translocate by a single base. Following Ambjornsson et al. [44], we find the change in the position of the centre of mass of the single stranded D N A coil located in the trans chamber upon a change in no of ± 1 . The relaxation time is then equal to the time required for the single stranded D N A coil to diffuse this distance. When fully dehybridized, the strands we consider are 76 bases long and have relaxation times ~ 10 s, _9 considerably smaller than the ~ lCT s typically required to translocate the D N A strand 6 by a single base. 6.3 Results and Discussion Equation 6.1 is used to compute the thermally averaged time for the D N A hairpin to first reach either set of absorbing states (n = do this, the (n = — NA, n{) and the (n = 2N + Afi 0 0 completely absorbing by setting P(-N , or (2N + —NA,UI) 0 oop + N ,ni L + A ^ n i = N). To = N) states are rendered ni)(t) = 0 and P(2N + A^ A loop + N , N) = 0 for L all t. Equation 6.1 is then dP(w, v, t) dt = /~^L{i,j\w,v)P(i,j,t) - P(w,v,t)^2L(w,v\i,j). ^(2N+N +N ,N) loop L (6.12) Defining A(t\(n' ,n[)) to be the probability that the system reaches an absorbing state Q between time t and t + dt, given that it started in state (n^n^) yields d_ M t \ ( « ) ) = E "dt yt(2N+N +N ,N) loop = L L(2N + N + N -1, P(2N + loop +N loop N\2N + N L loop + N ,N)x L N -l,N,t)+ L N J2 L{-{N A ni=0 ~ l)M - N ,m)P(-(N A A - l),n t). u (6.13) Chapter 6. First Passage Times of Driven DNA Hairpin Unzipping Name Sequence MFOLD Defect Free Defect-Bearing GCTCTGTTGCTCTCTCGCAACAGAGC GCTCTGTTGCTCTCTCGCAACTGAGC -31.9/c r B -27k T B 87 2-3 -26.0k T -24M T B B Table 6.3: D N A hairpin sequences [52] along with the hybridization energies predicted by M F O L D (T = 288) and the 2-3 free energy model. For more detail see [4]. Bases in the hairpin hybridization region are in bold, the defect bases are underlined, and the other bases form the loop segment of the hairpin. With the first passage time distribution A(t\(n' ,n[)) computed from the solution to Eq. 0 6.12, the mean first passage time to an extracted state is /•oo e(ni,n' ) 2 = / Jo tA(t\(n' n' ))dt 1} (6.14) 2 Finally, to find the thermally averaged mean first passage time, r, the results of Eq. 6.14 are summed over the initial states of the DNA, giving T y -pAG(o,o\i,j) e where AG(0, 0|i, j) is the free energy of state • (- ) b ib relative to state (0,0). In the experi- ments of Mathe et a/., the D N A hairpin was lodged and then held in the vestibule region for 0.5ms with a holding voltage of 20mV prior to the application of the driving voltage V [52]. So, when computing the value G(0,0\i,j) found in the Boltzmann weights in Eq. 6.15, we set the voltage equal to 20mV and sum over states where at least one hybridized base is found in the vestibule region. A simplex minimization code [3] was used to perform a least squares fit of Eq. 6.15 to the experimentally measured mean first passage times of the strands in Table 6.3 [52] (see Fig. 6.3(a)). While multiple local minima were found, all fits produced the same qualitative trends, and the best-fit parameters are given in Table 6.4. Fitting the model to both the defect-free and defect-bearing strands simultaneously shows that the 2-3 free energy surface provides a better fit than the surface produced by M F O L D . Specifically, when using the M F O L D free energies the fitted mean first passage times for the defectfree strand were found to exceed the experimental times, while the mean first passage times for the defect-bearing strand were less than the experimental times. As the primary Chapter 6. First Passage Times of Driven DNA Hairpin Energy Surface Ms- ) MFOLD 2-3 1.4 x 10 7.8 x 10 1 6 5 1.5 x 10 5.3 x 10 13 10 Q 0.39e 0.30e Unzipping 88 Sy 0.10k B omk B Table 6.4: Parameters from fitting to the experimental results for the mean escape time of D N A hairpins through an a-hemolysin pore of Mathe et al. [52]. The fit parameters sets were used to generate the mean escape time - voltage curves shown in Fig. 6.3(a). determinant of the mean first passage times is the hybridization free energies, these results indicate that in the context of this model, MFOLD's prediction of a 4.9k T B difference in hybridization free energy between the two strands may be an over-estimate. The value of the vestibule confinement entropy (s ) is also very small, suggesting that the D N A in v the vestibule region has considerably more conformational flexibility than D N A stored in the transmembrane region of the channel. We also note that the values for the charge per D N A base generated by the fits are larger than the values derived by Mathe et al. from single barrier fits by a factor of approximately 3. Nonetheless, these values for the partial charge are within the range found in the literature ( « 0.1 — 0.4e) [52, 53], for double stranded D N A dehybridization experiments of this type. Irrespective of which free energy surface was used, the curves in Fig. 6.3(a) predict an increase and peak in the mean escape time with increasing voltage at low voltages (below approximately 25mV). Figure 6.3(b) shows the probability (P erse) reV that the D N A escapes into the cis chamber and indicates that at low voltages the D N A will escape into the cis chamber almost exclusively. As the hairpin was initially inserted into the pore from the cis side of the membrane, we refer to escape into the cis chamber as reverse escape. In the high voltage regime, the hairpin dehybridizes and translocates into the trans chamber. To understand the reason for sharpness of the transition between the two escape mechanisms, we simplify the process of D N A escape and assume that the D N A strand faces a single kinetic barrier for each escape path. To determine the barrier to reverse sliding at low voltages in this simplified model, we assume the D N A does not undergo significant dehybridization prior to escaping into the cis chamber. Under these assumptions the maximal free energy is near state (-(N A - N ),0). L For forward escape, we set the free energy barrier equal to the hybridization free energy of the hairpin in the vestibule. Using Chapter 6. First Passage Times of Driven DNA Hairpin Unzipping 89 these assumptions to compute the barriers to reverse and forward sliding respectively, the probability that the D N A escapes into the cis chamber is approximately -(3AG(i,j\-(N -N )fl) k p ± A k e-0AG(i,j\-{N -N )fi) reverse r A L L -pAG(iJ\Q,N) + kf6 (1 + ^L (-P&G(-(N -N ),O\0,M))yl = where e A e L ^ is an arbitrary initial state, k is the effective attempt frequency for reverse r escape, and kf is the effective attempt frequency for forward escape. E q u a t i o n 6.16 has a sharp sigmoidal shape at the experimental temperature (288K). Solving E q . 6.16 for the voltage where P erse = 1/2 yields reV G + (2N - N )s T v V l / 2 v - \k T In {^^f) B ~ - A * T In ( £ ) \o\(N -N ) A ' L ( 6 - 1 ? ) where G is the hybridization free energy of the hairpin. A s reverse escape involves only sliding transitions, and the barrier to forward escape involves only dehybridization transitions, we set k = Li, kf = u>. In the case where the entropic penalty (s ) for storing r v bases i n the vestibule region is <C k , the numerator of E q . 6.17 is approximately equal B to G — k Tln(kf/k ). B r T h e difference i n the effective attempt frequencies then acts as a correction to the hybridization free energy. If the ratio kf/k becomes comparable to e 130 r then Vi/2 is approximately zero, and there is no transition between reverse and forward dominated escape. E q u a t i o n 6.17 suggests that the cross-over between reverse escape and forward escape will occur when the increase i n the electrical potential upon reverse escape becomes comparable to the hybridization free energy of the hairpin. Using the M F O L D free energy surface best fit parameters for the simultaneous fits (see Table 6.4), E q . 6.17 predicts transition voltages of 2 2 m V for the defect free strand, and 1 4 m V for the defective strand. These values are i n reasonable agreement with the predictions of the full stochastic model, which predicts 1 9 m V for the defect free strand and l l m V for the defect-bearing strand. W h e n the 2-3 free energy surface is used i n E q . 6.17 the predicted transition voltages for the defect-free and defective strand are 33mV and 29mV, while the full stochastic model gives 2 3 m V and 1 8 m V respectively. Note that E q . 6.17 used w i t h the 2-3 free energy surface gives higher transition voltages than Chapter 6. First Passage Times of Driven DNA Hairpin Unzipping 90 IOOOOF (a) 1000 k- 0.2 0.3 Driving Voltage (V) (b) i — Defect Free (MFOLD) — - Defect Bearing (MFOLD) •-• Defect Free (2-3) — Defect Bearing (2-3) 0.8 | 0.6 iC 0.4 0.2 0. 0 o.( (C) 2 0.01 0.02 Drr •ing Voltage (V) I • 0.03 r 0.04 • V=0.05 • V=0.15 • V=0.30 0.4 Figure 6.3: (a) The best fit predictions of the stochastic model when both the defect free and defect-bearing strands are fitted simultaneously. At voltages below the peak voltage, the D N A strands escape the pore primarily by sliding back into the cis chamber. At higher voltages, the hairpin completely dehybridizes and slides through the channel into the trans chamber. The inset shows that the voltage corresponding to the mean first passage time peak can be controlled by altering the length of the hairpin's poly-A tail, (b) The probability of the D N A strands escaping into the cis chamber, (c) The probability of having n bases pulled into the transmembrane region at the time the hairpin fully dehybridizes ( P a s t ) - Solid lines are the results for the defect-free 2-3 free energy surface, dashed lines are for the defect-bearing 2-3 free energy surface. Q Chapter 6. First Passage Times of Driven DNA Hairpin Unzipping 91 when used with the M F O L D free energy. This occurs despite the fact that the M F O L D hybridization free energy is larger than the hybridization free energies given by the 2-3 free energy surface, and is a consequence of the reduced ratio of u/p in the best fit parameters for the 2-3 free energy surface. The reduced ratio slows the rate of dehybridization relative to the rate of sliding, making reverse escape more probable. In the experiments of Mathe et al. mean escape times were measured at a minimum voltage of 30mV, and Mathe et al. did not find a peak in the escape time curve associated with the transition between reverse and forward dominated D N A escape. In the experiments of Nakane et al, a large Avidin anchor protein was attached to the end of the poly-A tail of the D N A strand, preventing reverse escape. This is equivalent to making the n = — N 0 A states in our model reflecting rather than absorbing. Making this change forces the D N A to dehybridize and escape into the trans chamber, producing an exponential rise in the mean first passage time at low voltages and eliminating the peak in the mean first passage time versus voltage curve. However, a switch between forward and reverse dominated escape is consistent with the measurement of the so-called entropic repulsion force by Nakane et al [53], who found that a minimum voltage of approximately lOmV was required to hold a double stranded D N A segment with a poly-A tail in a a-hemolysin pore. Similarly, Hendrickson et al. [37] found that to hold a completely single stranded D N A segment in an a-hemolysin pore, an applied voltage of at least 65mV was required. If we then assume that the reverse sliding of the D N A hairpin at zero applied voltage occurs on a time scale similar to the time for escape of single stranded D N A at zero applied voltage (a few milliseconds [41]), a transition between forward and reverse escape at a voltage in the range 11 — 23mV seems plausible. For a given hairpin, the location of the transition from the reverse to forward sliding regimes can be controlled by changing the length of the poly-A tail attached to the hairpin (inset Fig. 6.3(a)). Shortening the tail reduces the barrier to reverse sliding due to the applied voltage. As the poly-A tail is reduced, the voltage of the escape transition will then increase as suggested by Eq. 6.17. For voltages greater than the voltage associated with the maximum mean first passage time, the results of the stochastic model show a log-linear relationship between r and V . Chapter 6. F i r s t Passage T i m e s of D r i v e n D N A H a i r p i n Unzipping 92 In this model only the translation rate is influenced by the applied voltage. The reduction in the mean first passage time with increasing voltage is the result of inhibition of reverse sliding and inhibition of rehybridization due to the steric repulsion between the D N A strand and the membrane. This can be seen in Fig. 6.3(c) where -Pi t> the probability as distribution for having no hairpin bases pulled into the transmembrane region at the time the hairpin completely dehybridizes, is shown. As shown in Fig. 6.3(a), the discrete stochastic model also predicts that the mean first passage time will become relatively insensitive to changes in the applied voltage when V > 300mV. However, the applicability of the stochastic model at these voltages is limited by the barrier hopping model in Eq. 6.2. Specifically for a fixed value of ni, at large voltages Eq. 6.10 is approximately equal to A G r e When n 0 < 2N + N\ oop so that the transmembrane region is always filled with D N A , sliding of of the D N A is then described by the following master equation: = ~ p ( e - ^ v + l)P(n ,t) 0 + p e - ^ (6.18) where a is the distance moved by the D N A strand in a single sliding event, and A is the charge density along the D N A strand. In the limit of large applied voltage where exp(—/3|A|aV) -C 1, Eq. 6.18 becomes = p(P(n -l,t)-P(n ,t)). Q 0 (6.19) Equation 6.19, describes a system where only sliding toward the t r a n s chamber occurs, and where the rate of sliding is essentially insensitive to changes in the applied voltage. In our model the applied voltage does not directly influence the rate of dehybridization of the DNA hairpin. Thus, if the sliding of DNA becomes insensitive to the voltage as described in Eq. 6.19, the overall first passage time will also become insensitive to the voltage as seen in Fig. 6.3(a). Therefore, we expect that our results may not accurately model the influence of the applied voltage on the mean escape time at the highest voltages. We now consider the effects of the non-homologous base on the mean escape time. Chapter 6. First Passage Times of Driven DNA Hairpin Unzipping 93 At low voltages where the D N A escapes in reverse, Fig. 6.3(a) shows that there is no significant difference between the escape of the defect-bearing and defect free strands. This is consistent with the reverse escape mechanism where the D N A strands escape into the cis chamber and dehybridization is not the rate-limiting step in the escape process. Conversely, near the peak in the first passage times, the difference in the escape times between the defective and defect free strands is a maximum. Using the M F O L D free energy surface the maximum ratio T/Vdefect is approximately e , while the 2-3 free energy 42 surface gives a ratio of e . 18 In both cases the ratio is slightly smaller than the first passage time ratios that would be produced by a single barrier model of the translocation process ( e 49 and e , respectively). 2 We ascribe this difference to the coupling between the dehybridization and the D N A sliding into the transmembrane region that begins to occur near the peak voltage. Specifically, as the voltage is increased, reverse sliding of the strand is inhibited and the hybridization region is pulled deeper into the transmembrane region (See Fig. 6.3(c)). This has the effect of blocking the rehybridization of these bases, and Fig. 6.3(a) suggests this reduces the gap in mean first passage time between the defect-free and defect-bearing strands. Finally we look at the effects of varying the fixed parameters listed in Table 6.1 on our results. We begin by considering the effect of changing the hairpin length N. Changing N will change the total binding free energy of the hairpin G. Increasing N will then typically have the effect of increasing the mean first passage time, and increasing the voltage where the transition from reverse to forward escape occurs. Decreasing N sufficiently can eliminate this transition entirely. Equation 6 . 1 7 suggests this will occur when the total hybridization free energy is comparable to ksTlnikf/k ). r Using kf = u and k = p and r parameters from the fit to the 2-3 free energy, this suggests that the transition between forward and reverse escape would disappear when G « ll/c^T. If we assume an average hybridization free energy per base of 2.5/CBT, then this would mean the transition would disappear when N is approximately 4 or 5. At this D N A length, the stochastic model suggests the probability of escaping in reverse at zero voltage would approximately equal the probability of escaping in the forward direction. As was shown in the inset of Fig 6.3(a), scaling N A changes the location of the peak in Chapter 6. First Passage Times of Driven DNA Hairpin the mean first passage time curve. Increasing the value of N A Unzipping 94 means that on average more bases will be located i n the high voltage trans chamber i n the initial state of the hairpin. A s a result, there w i l l typically be a larger barrier to reverse sliding, the probability of escaping i n reverse at any given voltage will be reduced, and the voltage of the peak in the mean escape time will be reduced. We have tested for the existence of a peak in the mean escape time curve, and a transition from reverse to forward escape, for values of N A up to 60 using the defect-free 2-3 free energy surface and attempt frequencies. Varying has an effect similar to varying N . A W h e n Ni is reduced the number of D N A bases fully inserted into the trans chamber i n the initial state of the D N A typically increases. T h i s results i n a larger barrier to reverse sliding, and a reduction i n the voltage of the transition to forward escape. However using the defect-free 2-3 free energy, the transition remains for values of Ni as small as 2. Increasing A ^ , the number of single stranded bases that can fit i n the vestibule region, results i n a larger barrier to reverse sliding as the entropic penalty for reverse sliding increases w i t h N . v N v However because the fitted values of s v were relatively small, scaling has a minor effect on our results. The qualitative behaviour of the mean first passage time curve remained the same for values of to ksT, increasing N v up to 30. However when s is comparable v inhibits reverse sliding due to the increased entropic cost of filling the vestibule region. Altering A f loop , t h e number of D N A bases in the loop part of the hairpin, had very little effect on the mean first passage time results or the results for V i / . A s this quantity 2 only enters into the stochastic model when the D N A is fully dehybridized and nearly pulled through to the trans side, the results are relatively insensitive to variations in this quantity. Setting A f i o o p = 10 gives a mean first passage time curves and V\/2 values nearly indistinguishable from A ^ p = 5. Finally, i n order to reduce the number of fitting parameters i n our model, we have set UR, the rate of rebinding of a completely dehybridized hairpin, to zero. The simplest alternative to this assumption is to enforce detailed balance and set u R = u. This has the effect of dramatically increasing the m a x i m u m mean first passage times by more than 3 orders of magnitude compared to the results displayed i n F i g . 6.3, and increasing Chapter 6. First Passage Times of Driven DNA Hairpin Unzipping 95 the voltage of the transition between forward and reverse escape by more than 15 mV. Conversely, setting u R to a value much smaller than the sliding rate p has very small effects on the escape time curve displayed in Fig. 6.3 and produces increases in V~i/ 2 of approximately a millivolt. As UJ becomes comparable to p there is an increase of R a few millivolts in V i / , and an increase of about half an order of magnitude in the 2 maximum mean escape time. However, the predicted mean first passage times at voltages greater than 30mV where the majority of the experimental data was taken, do not change significantly for the 2-3 free energy surface. As we do not have any data outside this range, we have applied the assumption made in single barrier models of D N A escape previously used to study the D N A translocation process, and set u> = 0. This could easily be revised R in light of new experimental results. 6.4 Summary We have analyzed the processes of translocation of D N A hairpins through narrow transmembrane pores. Using a simple two dimensional stochastic model, we produced estimates for the mean first passage time which suggest that for voltages below approximately 25mV, the D N A hairpin can escape the transmembrane pore without completely dehybridizing and without traversing the membrane channel. When escaping the transmembrane pore in the reverse direction the influence of the D N A sequence on the escape time is minimal. At voltages greater than approximately 25mV the D N A strand completely dehybridizes and translocates into the trans chamber. The switch between escape into the cis and trans chamber results in a peak in the mean escape time curve. The peak, and the transition between escape into the cis and trans chambers, occurs when the free energy to dehybridize the hairpin approximately equals the free energy penalty for pulling the poly-A tail from the high voltage trans chamber to the low voltage cis chamber. Additionally, near the peak in the mean first passage time curves the two D N A strands examined exhibited the largest differences in mean escape time. This suggests that in a device that used mean escape times to distinguish between strands of different sequences, it would be desirable to adjust the driving voltage to coincide with the peak of the mean first passage Chapter 6. First Passage Times of Driven DNA Hairpin Unzipping 96 time curve of at least one of the D N A strands in the sample under examination. Between approximately 25mV and 200mV, the mean escape time had the approximately log-linear dependence observed in experimental measurements of the mean translocation time. 97 Chapter 7 Water Adsorption in Ion bearing Nanopores 7.1 Background The behaviour of water and electrolytes in confined geometries is a topic of significant current interest with applications to many areas of bio- and nanotechnology. In particular, the interaction of water and ions in the highly restrictive confines of cellular transmembrane channels [54-56, 59, 61] is of extreme importance to the functioning of biological systems. Of tremendous interest in the biophysical context, is the relationship between ion solvation and the dramatic selectivity for ions of different species displayed by biological channels [92]. Similarly, the interaction of water and molecular vessels such as carbon nanotubes [93] and fullerenes is receiving significant attention in an aim to develop nanoscale fluid conduits, and storage devices [61-64]. Most existing work on the behaviour of both pure water and electrolytes in nanoscale confined geometries, considers confinement in infinite, slit-shaped or cylindrical pores. In both cases, computer simulations have been extensively employed to investigate the water and ion behaviours. We will focus on cylindrical pores where an early molecular dynamics study of Lynden-Bell and Raisiah [60] considered the behaviour of single ions solvated with explicit water molecules in narrow tubes with strictly repulsive walls. The results of these simulations show significant water structuring in the tube at room temperatures (298 K) induced both by the restricted geometry and the ion. Additionally, equilibrium ion mobilities were computed and were found to be 0.25 to 0.5 times those of the ion solvated in the bulk. Similar studies have investigated properties such as non-equilibrium ion mobility [94], ion selectivity [57], and ion pairing in narrow pores [95]. However, in all these studies the equilibrium density of water in the pore was never computed, and instead was set equal to the bulk value. A natural means of calculating this density is Chapter 7. Water Adsorption in Ion bearing Nanopores 98 through Monte Carlo simulations in the grand canonical ensemble where water molecules are free to enter and exit the pore. While very recent work by Striolo and collaborators, [96-99] and Liu and Monson [100], employ simulations in the grand canonical ensemble to examine the adsorption of water into the interiors of carbon nanostructures, the pores in these simulations do not contain ions. In this chapter, the solvation of ions bound in nanometer-scale cylindrical pores will be examined using grand canonical Monte Carlo simulations. Following previous work, we include only one ion in the pore at a time [57, 60], resulting in a violation of charge neutrality in the pore interior. Charge neutrality violations have been observed in simulations of electrolytes in highly confined nanopores [101-103], where wall charges are present. Additionally it has been argued [95] and observed [58] that neutrality violation is also possible in narrow pores in the absence of wall charge. While local neutrality in the pore interior can be violated, global charge neutrality is always preserved. Thus the simulations in this chapter can be interpreted by considering the ion to be part of the pore itself, and where the counter-ion is trapped in the water bath. The results of this chapter focus on the adsorption isotherms at room temperature and determine the effects of the ions on the filling process. We then consider the water structure in the channel at equilibrium at various values of the bulk chemical potential, and based on an examination of the water structure, discuss the mechanism driving adsorption. 7.2 7.2.1 M o d e l and Simulation Methodology Molecular Models The implementation of a molecular simulation of water and ions in a confining cylindrical pore requires three basic components; models for the water molecules, ions, and the pore, a Hamiltonian specifying the interactions between these components, and a mechanism for evolving the state of the system in a manner consistent with the Hamiltonian. We begin by considering the water molecules which we model using the well-established Simple Point-Charge/Extended (SPC/E) interacting site model for water [60, 96, 106, 107]. The S P C / E model for a water molecule can be described as a rigid body of three Chapter 7. Water Adsorption Species cr(nm) e(kJ/mol) SPC/E C Na K Ca c r 0.3169 0.3283 0.2876 0.3250 0.3019 0.3785 0.3143 0.6502 0.3891 0.5216 0.5216 0.5216 0.5216 0.6998 + + 2 + in Ion bearing q 0 Nanopores 9(e) = -0.8476, q 0 1.0 1.0 2.0 -1.0 -1.0 H 99 = 0.4238 Table 7.1: T h e potential parameters (at T = 298K) used for the different species [60, 96, 104]. A l l Lennard-Jones parameters are for interactions between the listed species and S P C / E water. Cross-interaction parameters were computed using Lorentz-Berthelot mixing rules [105]. T h e carbon parameters were used for the wall potential given in E q . 7.2. point charges; a negative point charge corresponding to the water oxygen and two positive point charges corresponding to the water hydrogens. T h e water molecule interacts with the surrounding environment v i a the Coulombic interactions associated w i t h the point charges, and through a Lennard-Jones potential [105] centred on the oxygen atom. B o t h the partial charges and the Lennard-Jones parameters were originally determined by fitting the results of isobaric molecular dynamics simulations to experimental values for the bulk density and dielectric constant of liquid water at 298K [106]. T h e values for these parameters are listed i n Table 7.1. T h e model used for the ions consists of a point charge that interacts v i a a Coulomb potential and a concentric Lennard-Jones interaction. T h e parameters for the various ionic species used i n this study are also given i n Table 7.1. T h e geometry of the confining cylindrical surface is displayed i n F i g . 7.1(a) where the symmetry axis of the cylinder is taken to be the z-axis. T h e cylinder is modeled as featureless smooth wall using one of two potentials. Following L y n d e n - B e l l [60] we employ a purely repulsive exponential potential 0waii = 0oe- ; ( K - , r ) (7.1) to model strictly hydrophobic tubes. In E q . 7.1, R is the radius of the confining cylinder, and r = \ J X 2 + V 2 LS the radial position of the particle. Note that as there is no dependence on z in E q . 7.1, the confining cylinder is homogeneous i n the z direction. T h e value of Chapter 7. Water Adsorption in Ion bearing Nanopores 100 / in Eq. 7.1 sets the effective length of the wall interaction and here we set / = 2.0nm _1 resulting in a very short-range potential that provides a good approximation to a perfectly hard wall. While Eq. 7.1 is simple to implement and compute, it does not provide a particularly realistic model of the interior of a cylindrical nanostructure such as a carbon nanotube. To enhance the realism of our simulations, we also consider the case where the walls of the cylindrical pore interact with the confined particles using a smooth potential derived from the well-known Lennard-Jones interaction [105]. To compute this potential, we consider an infinitely long cylinder with a constant and continuous surface density of Lennard-Jones centres, and then integrate over the surface of the cylinder. This gives [108] OO /»7T / O / w a l ? / = eno 2 2 7r OO J — 7T 63 F 32' \ -9/2,-9/2,1, 3F(-3/2,-3/2,l,(^) ) 2 In Eq. 7.2, n is the surface density of Lennard-Jones centres on the bounding cylinder, F is the standard hypergeometric function, and TZ is the vector between the particle at f and a point on the cylinder wall. Note that while n is determined entirely by the structural details of the bounding cylinder, it ultimately acts only as a correction to e, the LennardJones interaction strength. Thus, while the strength of the potential in Eq. 7.2 is strongly dependent on the nature of the bounding cylinder we choose, the qualitative characteristics of the bounding potential are not. In this chapter we set n — 4 0 n m , the approximate -2 atomic surface density of single-walled carbon nanotubes [93]. However, given that the qualitative characteristics of Eq. 7.2 do not depend on n, we expect our conclusions will hold qualitatively for a variety of cylindrical pores and channels. Plots of Eq. 7.2 for several cylinder radii can be viewed in Fig. 7.1(b). The Lennard-Jones parameters (o and e) used in Eq. 7.2, are derived from the parameters given in Table 7.1 by applying the Chapter 7. Water Adsorption in Ion bearing Nanopores 101 Lorentz-Berthelot mixing rules [105]. Finally, in Fig. 7.1 and Figs. 7.6 through 7.13, the bounding cylinder is displayed as atomistic rather than smooth. The bounding cylinder displayed in these figures should be treated as an aid to understanding the qualitative features of the structures of the adsorbed particles, and not a literal representation of the source of the bounding potential given by Eq. 7.2. Having defined the molecular models used for all parts of the system, we can write the Hamiltonian, H — f/ kin + <7 conf . In the Hamiltonian, the kinetic term is derived from the standard expressions [109] for the kinetic energy of a point particle representing the ion, and for a rigid body representing an S P C / E water molecule. The orientation of each water molecule in space is described by a set of Euler angles [109]. The total kinetic contribution to the Hamiltonian is then [110] nw Uconf = 0wall(7"o\ &WWI, ^WWl) + ^ w a l l t j j , VlWh ^IWl) + i=l n w n / w i=l j=i+i \ // \\'o >o\j sr^ sr^ Qm \\v sr~^ 2 E «=i E - me{0,H H } lt 2 m 47re |r-/-r. 0 r q n u w 6 o I o ne{0,H H } I' i=l j=i+lme{0,Hi,H } n \ 6 2 m ' n I Chapter 7. W a t e r A d s o r p t i o n i n Ion b e a r i n g Z = — -- R R R R = = = = 0.44nm 0.54nm 0.64nm 0.74nm Nanopores 102 0 T—i—i—nn—i—i—I—i—i—i—i—i—m—1~ / I < I I I \ I o.i _L J_ 0.2 0.3 0.4 r (nm) 0.5 Figure 7.1: (a) The geometry of the confining pore. Unless otherwise noted L — 1.5nm. (b) The wall potential as a function of radial position. The radius at which the water-wall potential energy is equal to k T is taken to equal the physical radius of the pore. B Chapter 7. Water Adsorption in Ion bearing Nanopores In Eqs. 7.3 and 7.4, quantities with the subscripts O, Hi, H 2 103 refer to the water oxy- gen, and hydrogens respectively, while quantities with the subscript I refer to the ion. I\, 1%, and I3 are the principal moments of inertia of an S P C / E water molecule about the molecule's centre of mass, and pj is the centre of mass momentum of particle j. (to, 9, <p) are the Euler angles used to describe the orientation of the S P C / E water molecules about their centres of mass, and P^e^} are the conjugate momenta associated with the Euler angles of a water molecule i. The quantities (oiwi, tiwi), (owwi £ww)y {owwi, wwi), e (OIW,£IW), and are the Lennard-Jones parameters for the ion-wall, water-wall, ion-water, and water-water interactions respectively. With the Hamiltonian defined, the last remaining parts of the model to be specified are the boundary conditions applied along the axial dimension of the confining cylinder. One possible choice is to apply periodic boundary conditions along the symmetry axis of the pore, while simultaneously applying the minimum image convention [105] to the interactions between the ions and water molecules in the pore. Under this convention, a molecule or ion only interacts with the nearest periodic images of the other particles in the system. As a result, the Lennard-Jones and electrostatic interactions specified in Eq. 7.4 are effectively cut-off for distances longer than L. Note however, that the interactions between the wall and the confined particles are not truncated when using Eq. 7.2. In principle, by applying periodic boundary conditions to the ion-bearing pores we are not simulating pores with a single ion, but infinite pores with a finite ion density. If all the interactions in the system were computed, this would result in the simulation cell having an infinite energy. However, with the truncation on the inter-particle interactions from the minimum image convention, the ion-ion interactions are reduced to zero and the energy of the simulation cell remains finite. The use of periodic boundaries with the minimum image convention, does however lead to questions about the physical interpretation of the simulation results. A simple interpretation of the simulations with the minimum image convention and periodic boundary conditions is that they approximate a finite pore of length 2L. The drawback to this interpretation can be seen if we consider that the ion will induce a polarization in the water molecules with values of the z coordinate greater than the ion's, and the opposite polarization in water molecules with values of z less than Chapter 7. Water Adsorption in Ion bearing Nanopores 104 the ion's. Applying periodic boundary conditions will then result in a spurious interface between water molecules with opposing dipole orientations. To avoid this induced polarization interface, we could create a real interface by capping the ends of the nanopore and simulating a closed system. However, this naturally leads to questions about how water could enter the pore. Additionally, we would like to extrapolate these results to finite open systems, such as biological ion channels, where all the water at the ends of the channel would be free to interact with water in the bulk. In this case, a capped pore that eliminates the water-water interactions at the pore ends would arguably be no better an approximation than the pore under periodic boundary conditions. Given all these factors, we elect to follow previous work on ion solvation in pores and in the bulk [60, 104] and apply periodic boundary conditions with the minimum image convention. The effects of this choice on our results will be explored later in the chapter. 7.2.2 Simulation Methodology Similar to the case of the T A S E P discussed in Chapters 2 through 5, where the steadystate lattice probabilities were sought, here we seek the equilibrium probabilities of finding the water-ion-pore system in its various possible occupancy states. Given the Hamiltonian, a state of the water-ion-pore system is defined by the number of water molecules in the channel (nw), the position of the centre of mass of each of the water molecules (ft), the orientation angles of the water molecules Qi = (uJi,6i, (pi), and the position of the ion (fi). In the long-time limit, a system described by the Hamiltonian given in the previous section, and coupled to a heat bath at temperature T and a water reservoir with chemical potential pw, will come to thermal equilibrium. The state probabilities for the set of unlabelled particles will then be given by the Boltzmann distribution. The grand canonical partition function of the system is given by Chapter 7. Water Adsorption poo 0 0 Z G = in Ion bearing V (nwl)- ^)-^h~ 1 / 3 „n w = 0n ^-00 n d dn 3 3 Pl p2n pir .../ JO p2rr / / Jo Jo Nanopores pL .../ J-L p2-n 105 pR / / Jo Jo -MH-,* n ) e w w x w Y[ d^dpMdp^dnid ^ 3 i=l (7.5) After integrating out the momentum terms, the partition function can be written (2A e e e- ) Jo nw nw=Q n\v=0 n w 1 2 i JO J-LJO JO JO w Y\ sm(8i)dQid r d r , : i 3 i I 1=1 . where U { is the configurational potential energy of the system, j3 = 1/ksT, con (7.6) A / is the thermal wavelength of the ion, Aw is the translational thermal wavelength of the S P C / E water molecule, and ( © 1 , 6 2 , ©3) are the rotational thermal angles of the S P C / E water molecule ^/2iimikBT 1 ) A = ©* = h \j2irniwksT' h , yJ2ixIik T B In the expressions for the thermal wavelength and thermal angles, mi is the ion mass, m-w is the water mass, and f is the i th principal moment of inertia of the S P C / E water molecule. The equilibrium probability of any particular state of the system is then Chapter Pe (n , q w 7. Water Adsorption {ji}, {f2/},f/) in Ion bearing = ZQ A7 (e-^2A ^eie e )3 1 3 2 3 n M / Nanopores e- / 3 t / p° t ( r i M / 106 ' ' ^ ^ x { r l } { } ) J^[(i r sin(6 )(ir2 (i r/. i=l 3 l i 3 i i (7.7) While Eq. 7.7 in principle provides all the information needed to compute any equilibrium property of the system, directly computing the state probabilities is not feasible due the enormous computational cost of computing the high-dimensional integrals in the partition function. This is directly analogous to the situation described in Chapter 2 where the large state space of the T A S E P made direct solution of the occupancy equations impractical. As was the case with the TASEP, we take an indirect route to compute the probabilities of Eq. 7.7 by defining a stochastic process that has a long-time probability distribution equal to Eq. 7.7. We then use a Monte Carlo algorithm to evolve the system toward this distribution. While the confined electrolyte considered in this study can take on any of an infinite and uncountable number of states, for the sake of notational clarity we will consider a simplification where the number of states is treated as countable. Other than making the notation less cumbersome, this simplification has no other effect on the description of the simulation methodology [111]. With this modification, the evolution of the stochastic process used to model the confined electrolyte is then governed by a discrete-state Master equation P n+1 = MP , (7.8) n where P is a vector of state probabilities, and M is a matrix of transition probabilities between the various states of the water-ion-pore system. The subscript n refers to the number of iterations of the stochastic process that have been completed. In writing Eq. 7.8 we have implicitly enumerated all the possible states of the confined water-ion system and labelled the states with a single index. Element of M is then the probability that the system in state j transitions to state i and is a product of two quantities, M j = t a^Aij. Here aij is the probability that a transition from j to i is attempted, and Aij is the Chapter 7. Water Adsorption in Ion bearing probability that the transition is successfully completed. Nanopores 107 Regardless of the stochastic dynamics we choose, the transition probabilities must satisfy four conditions. First, they must conserve probability which requires that X X = 1 - ( 7 - 9 ) i They must also satisfy the condition of detailed balance Mii = M^l (7 m) Additionally, the stochastic dynamics must also be formally regular, meaning that given an arbitrary state i, there must be a finite number of transitions that will move the system to another arbitrary state j. Finally, we require that M y be a function of (i, j) only, and thus be unaffected by state transitions previously attempted by the system. If all these conditions are satisfied, the transition probabilities define a so-called Markov Process that is guaranteed to evolve to the equilibrium probability distribution given by P e q [68]. We make use of the well-known Metropolis algorithm [66] to define the acceptance probabilities (Aj), for the following set of transitions 1. Water insertion: A water molecule is inserted into the cylinder by placing the centre of mass at a randomly selected location, and then rotating the molecule to a random orientation. T h e attempt probability is 5 sm(6)dxdydzdud6d(p » wv~ °V = where 8 is an arbitrary value in the range [0,0.5). T h e acceptance probability is then Aj = e^ m a x ( { c ^ 0 ) 1 where 47r VV 4TT 2 / Cij = (Ui - Uj) - 2 \\ - k T\n (7-12) . B \A- QiQ2®3(n (]) w w + i)y In E q . 7.12 (Ui — Uj) is the difference i n the total potential energy after and before Chapter 7. Water Adsorption in Ion bearing Nanopores 108 the insertion, V is the total internal volume of the channel, and nw(j) is the number of water molecules in state j. 2. Water deletion: A water molecule in the cylinder is selected at random and removed from the cylinder with attempt and acceptance probabilities given by Oij = S/n , (7.13) w and A tj = , dij = e -/? max ( K ; > ° } ) where , (Ui -Uj)+ m T r M . rp, f^w^2e n (j)\ - k TIn I 3 p w w B J. (7-14) 3. Water rotation: A water molecule currently inside the channel is selected at random, and is rotated through a randomly selected angle about a randomly selected axis. The attempt and acceptance probabilities are / Aij = 1-25 \ {sm(6)dujd6d(f)\ -P™^m-u )fi}). e (. ) 3 7 16 4. Water Translation: A water molecule currently inside the channel is selected at random and is translated a randomly selected distance in a randomly selected direction ( 1-25 \ (dxdydz\ 5. Ion Translation: The confined ion is translated a randomly selected distance in a randomly selected direction Chapter 7. Water Adsorption in Ion bearing Nanopores 109 (7.19) (7.20) The factors of one half in the expressions for in the water translation and water rotation moves are purely a result of the way the simulation selects between translations and rotations. In the simulation, a translation or rotation event of some kind occurs with probability 1 — 28. The chance of any particular particle being selected is (n + l ) - If - 1 w the selected particle is a water molecule, either an attempt at a linear translation is made with probability 1/2, or an attempt at a rotation is made with probability 1/2. If the selected particle is the ion, the probability of a linear translation attempt is 1. The attempt and acceptance probabilities defined above have been explicitly chosen to ensure that the detailed balance condition is met. The regularity of the stochastic process can be seen by noting that states can only differ by the number of water molecules, the position of the water molecules and ion, and the orientations of the water molecules. Arbitrary differences in riw can be eliminated by repeated insertions or deletions, while the rotational and translational transitions can eliminate arbitrary differences in configuration. Finally, note that if an attempted state transition is rejected, the state of the system does not change. Thus in all cases the diagonal elements of the transition matrix are given by M 3j = 1 — J2i,i^j ^ij> thereby ensuring that Eq. 7.9 is satisfied. In our study of the water-ion-pore system we set 8 = 0.25. Additionally, the range of the random translations and rotations permitted in a single transition attempt were adjusted so that these transitions were accepted approximately 50% of the time. In order to ensure the system had equilibriated, 10 transition attempts were run before 9 data collection began. The approach to equilibrium was monitored by considering the variation in the total energy and total pressure as the number of attempted transitions increased. Once any systematic variation in these quantities with respect to the number of attempted transitions was eliminated, and the quantities randomly fluctuated about well-defined mean values, the system was considered to have reached equilibrium. After Chapter 7. Water Adsorption in Ion bearing Nanopores 110 equilibrium was achieved, statistics were taken for as many as 3 x 10 transition attempts. 9 Unless otherwise noted, the period of the cylindrical pores was set to 3.0nm, a length which is consistent with the typical length of biological transmembrane channels [87]. The cylinders we examine have centre to wall, or geometric, radii of 0.44, 0.54, and 0.74nm. 0.64 However, due to the extended repulsive core of the carbon atoms in the cylinder wall, the physical radii of the pores are considerably smaller than the geometric radii. Setting the physical radii equal to the one k T B point of the water-wall interaction potential (Eq. 7.2) gives a physical radius of approximately 0.16nm for the pore with a geometric radius of 0.44nm. Each O.lnm increase in the geometric radius then produces a O.lnm increase in the physical radius. Note that all simulations in this chapter were carried out at a temperature of 298K. Additionally, the energies in this chapter are typically quoted in kJ/mol where, unless otherwise noted, this refers to the energy in kilojoules of a mole of pores in the given state. 7.3 7.3.1 Results and Discussion Adsorption Isotherms Figure 7.2 displays the adsorption and desorption isotherms for water in ion-free, weakly hydrophilic pores. The vertical line in the figure defines the chemical potential for saturated water vapour at T = 298A" (p w —58.5kJ/mol [100]). For the adsorption curves the simulations were initially started with no water in the pore. The initial condition for the desorption branches was an arbitrarily selected water configuration from the point on the adsorption curve with the largest value of p . w Consistent with earlier work, the pure- water isotherms exhibit a sharp filling transition for all pore sizes greater than 0.44nm. Note the value of pw for filling decreases when the radius is increased from 0.44nm to 0.54nm, and then increases for all larger radii. As noted in earlier studies [96], this nonmonotonic behaviour is the result of the mean water-wall interaction energy becoming less favourable as the pore radius increases, while the penalty for disrupting the preferred tetrahedral structure of liquid water decreases as the radius is increased. Also, in agreement with previous studies [96, 97, 100], the pure-water adsorption-desorption curves for Chapter -62 7. Water Adsorption -61 -60 -59 in Ion bearing -58 -57 -56 Nanopores -55 111 -54 Figure 7.2: Adsorption and desorption isotherms for water and non-ion-bearing pores. The solid lines correspond to the adsorption isotherms, while the dashed lines correspond to the desorption isotherms. Hysteresis was not observed in the 0.54 and 0.44nm pores. Lines are drawn as guides only. At the density of liquid water at 298K (l.Og/cc) the expected number of waters adsorbed in the 0.44nm pore is 8, 21 in the 0.54nm pore, 41 in the 0.64nm pore, and 67 in the 0.74nm pore. Chapter 7. Water Adsorption in Ion bearing Nanopores 112 the 0.64 and 0.74nm pores show significant hysteresis that increases with increasing pore radius. Figure 7.3 shows the isotherms for pores with geometries identical to those of Fig. 7.2, but now containing a mobile ion. The most striking effect of the presence of the ion is a reduction in the value of u, w at the filling point. This reduction is a consequence of the decrease in the chemical potential of the confined water due primarily to the favourable electrostatic interaction between the water and ion. Also note that even at low chemical potentials, the ion-bearing pores always retain a non-zero number of water molecules that form a loose cluster around the ion. This is consistent with previous studies of water nucleation near ions [112] that have shown that an ion's first hydration shell is retained as the ion is moved from liquid to vapour. Additionally when an ion is present in the pore, the sharpness of the filling as a function of u. is visibly reduced, particularly in larger w pores with the doubly charged calcium ion. This indicates a change in the filling process from a sudden collective effect in the ion-free case, to a more gradual growth mechanism in the presence of the ion. This change in the filling mechanism is also evidenced by the absence of significant hysteresis in the adsorption-desorption curves for the ion-bearing pores. Using the volume computed from the physical radius, at the bulk density of liquid water at 298K and latm pressure (~ 1.0g/cm ), the 0.44nm pore would have approx3 imately 8 water molecules, the 0.54nm pore would have approximately 21, the 0.64nm pore would have approximately 41, and the 0.74nm pore would contain approximately 67 water molecules. As can be seen in Figs 7.2 and 7.3, when filled both the ion-laden and the ion-free pores have a water density exceeding the bulk density of water at 298K and latm pressure. Note that in this context, a filled state refers to any state of the pore where the majority of the interior pore volume is occupied by water molecules. Another interesting feature of the filling isotherms of the ion bearing pores, is the variation in the number of water molecules in the filled pore with ionic species. Consider the isotherms for the monovalent cations N a and K . In general, pores carrying a sodium + + admit more water molecules in the filled state than pores carrying a potassium, and this effect is largest in pores of small radii. This is entirely consistent with the smaller van der Chapter 7. Water A d s o r p t i o n i n Ion b e a r i n g Nanopores 113 Figure 7.3: Adsorption isotherms for water into ion-laden pores. The horizontal line gives the expected number of water molecules in the physical volume of each pore at the ambient density of liquid water 1.0g/cm at 298K), while the vertical line gives the chemical potential of saturated water vapour at 298K ( « -58.5kJ/mol). 3 Chapter 7. Water Adsorption in Ion bearing Nanopores 114 Waals radius for sodium (see Table 7.1) that both reduces the excluded volume associated with the ion and allows water molecules to approach the ion more closely leading to a universally more favourable ion-water electrostatic interaction. However, comparing the monovalent cations to chloride shows a striking variation in the density at the filled state. In the narrowest 0.44nm pore, the chloride-bearing pore has a lower density than the cation-bearing pores consistent with chloride's monovalency and comparatively large van der Waals radius. However, in the larger pores, the presence of a C P (or F~) results in an increase in the filled density compared to the monovalent cations, and approximately matches the density induced by the presence of a C a . 2 + 7.3.2 Filled State Structures We begin by considering the structures formed in the pure water pores in the filled state (uw = — 56.0kJ/mol). Consistent with earlier results [60, 96] we see a single line of water molecules in the narrowest 0.44nm pore, and the formation of water layers in pores of larger diameter (Fig. 7.4). Particularly interesting is the water in the 0.54nm pore that forms into a relatively stable hydrogen-bonded structure. Assuming a hydrogen bond forms when TOH (the oxygen - hydrogen distance) is < 0.25nm, and the O H O angle is greater than 130° [60], each water in the 0.54nm pore participates in approximately 3.1 hydrogen bonds. The water structure in the 0.54nm pore is consistent with the empirical hydrogen bond rules suggested by Koga et al. [113], with the hydrogen bonds lying in a common cylinder plane sharing a common orientation (i.e., all clockwise or counter-clockwise oriented). The stacking pattern of unit cells proposed by Koga is not strictly observed, but with the energy difference between various stacking patterns on the order of 0.05rCBT [113], this is not unexpected. The 0.64nm and 0.74nm pores show a layering effect where water coats the inner surface of the pore. In the 0.64nm pore at p. w = —56.0kJ/mol, approximately' four water molecules are also found in the centre of the pore. In the 0.74nm pore a second water layer begins to form in the centre of the pore with a pore-centre-to-water-oxygen radius of approximately 0.15nm. Chapter 7. Water A d s o r p t i o n i n Ion b e a r i n g Nanopores 115 Figure 7.4: Characteristic water structures in pores with geometric diameters of (a) 0.44nm, (b) 0.54nm, (c) 0.64nm and (d) 0.74nm. The line of water molecules in the 0.44nm pore is nearly perfectly hydrogen bonded with an average of 1.8 hydrogen bonds per water molecule, pw = —56.0kJ/mol. Chapter 7. Water Adsorption in Ion bearing Nanopores 116 The water-oxygen-water-oxygen radial distribution functions in various pores are displayed in Fig. 7.5, where the radial distribution function is given by - r £ E i=l j = ljy£i 5 3 ( fi ~ i f _ ^) e-^ p o t ( " r / , i ? 1 - - ^ ) d r fj x d r sin e dQ d TZ. —\ • 3 / 3 3 fc k k k (7.21) As before, 77 in Eq. 7.21 refers to the position of the ion, and here fj is the position of the centre of mass of the i th water molecule. Figure 7.5 shows the addition of an ion to the pores at the same water chemical potential produces a significant disruption in the water structure for the narrower pores, and comparatively little in the pores with larger radii. Consistent with its larger charge, the addition of a C a 2 + ion results in larger deviations from the pure water structure than does the addition of a monovalent ion. As observed previously [60], the addition of an ion creates a long-range orientational ordering in the pore with the water dipoles preferentially aligning with the ion's field. A n example of this ordering can be seen in Fig. 7.6(a). As a consequence of this alignment, pores with periodic boundary conditions possess an interface at the period boundary between water molecules with dipoles oriented in opposing directions relative to the pore axis. As we will show later in this chapter, the opposing dipole orientations induced by the ion do not dramatically perturb the equilibrium water structure in the pore. 7.3.3 Adsorption In Ion-Bearing Pores Figures 7.6 through 7.13 show sample configurations of water in pores with K Cl _ at different values of pw along the filling curve. + and In each case the configurations displayed include the most probable number of water molecules at the given chemical potential. As evident in Fig. 7.3, at low chemical potentials the ion-bearing pores all contain approximately 4 to 10 water molecules strongly localized near the ion. As ionfree, but otherwise identical, pores are empty at chemical potentials below —59kJ/mol, the water cluster is stabilized in the pore by the favourable electrostatic interactions with Chapter 7. Water Adsorption in Ion bearing Nanopores 117 Figure 7.5: Water Oxygen - Water Oxygen radial distributions functions for pores containing pure water, and pores containing single ions (p = —56.0kJ/mol). The ions produce deviations from the pure water structure that are most pronounced in pores with smaller radii, and increase with the size of the ion's charge. In the 0.74nm pore, the perturbation of the water structure resulting from the ion is small. w Chapter (a) 7. W a t e r A d s o r p t i o n i n Ion b e a r i n g N a n o p o r e s (b) (c) (d) (e) 118 (f) Figure 7.6: Filling in the 0.44nm pore with a C P ion. (a) p w — —66kJ/mol (n = 4) (b) p w = —64kJ/mol (n = 6) (c) the low density state at p w = —62kJ/mol ( n w — 7) (d) the high density state at p = —62kJ/mol ( n w = 11) (e) p w = —59kJ/mol (n = 12) (f) p w = —56kJ/mol ( n w = 13). At p w = —56kJ/mol states with 12 water molecules similar to the state shown in (e) and states with 13 water molecules are equally probable. w w w w the ion and not the interaction with the wall potential. A bulk chemical potential of approximately —90kJ/mol corresponding to extremely rarefied water vapour at a density of ~ 1 0 m o l / m , will completely empty the pore of water. _6 As 3 the reservoir chemical potential is increased, Figs. 7.6 through 7.13 show that the water cluster around the ion grows but remains highly localized. Figures 7.14 and 7.15 show the water-ion, water-wall, water-water, and total configurational energy for the monovalent ions. When comparing ions with the same charge, smaller ions produce more favourable water-ion interaction energies, and less favourable water-water interactions. The wall interaction is relatively insensitive to ion size, but the ion-wall interaction does slightly favour the larger ions. Figures 7.16 and 7.17 show the radial probability densities for the various pores carrying a C P and a K + ion, respectively. Note that in all cases there is a strong tendency for Chapter (•••••I !•••••* '•••••Ii< '•••••HI !••••••> • ••••II •••PI' 7. W a t e r A d s o r p t i o n i n Ion b e a r i n g 88&1 •*'JMHj '•••••I '•••••I >•••••!! '•••••II ••£••»• M*J«I< • ••••• •I ' ••?*")• '•••*•» •rr*«i •••••i •••••• •••••li ••••••• '•••--••I •••••II •••••II •••••! •••••i •••••i •••••• : : : : : : : '•••••III ( • • • • • I I I • • • • • • I I I ilPi •••••II • •^••Il( • • • • • • I • •^••l • • • • • i i !•••••!• !•••••!• _ JfMli '••W««H '•••> ,)|l< '••••»•••( was, • ••"»••!« • • • • • • I fe >•••••!• • ••••II • ••••II • • • • • • I «•••••••< • • • • • i i (•••••II •••••II Nanopores 119 I M I I I II > • * • ( " • •••W •n»«i •-•••i .<ri Hl.<r| ••fcOM • • * H * •••Mi •^*o* •!•••! p ••'•••« ••••• !••*••» !•••*••!< ••••ii € ^ $?) (s) (§) (J) Figure 7.7: Filling in the 0.44nm pore with a K ion. (a) p = -66kJ/mol ( n w = 5) (b) p = -64.0kJ/mol (n = 7)(c) p = -62.0kJ/mol (n = 9) (d) P w = —61.0kJ/mol ( n w = 10) (e) p — —59.0kJ/mol (n = 13) (f) p w = -56kJ/mol (n = 14). + w w w w w w w w Chapter (a) 7. Water Adsorption (b) in Ion bearing Nanopores (c) (d) 120 (e) (f) Figure 7.8: Filling in the 0.54nm pore with a C P ion. (a) p w = -66kJ/mol ( n w = 6) (b) p = —64kJ/mol (n = 8) (c) the low density state at p w = —61.1kJ/mol ( n w = 15) (d) the high density state at p = —61.1kJ/mol ( n w = 31) (e) P w = —59kJ/mol ( n w — 34) (f) p = —56kJ/mol ( n w = 40). w w w w the water oxygens to layer near the minimum of the wall potential. In simulations using the wall potential of Eq. 7.1, layering is not observed prior to filling in pores with physical radii greater than 0.3nm; this indicates that geometric confinement alone is not sufficient to produce layering before the pore fills. After filling however, layering is observed in pores using Eq. 7.1, consistent with the observation of the formation of polygonal ice structures in nanotubes driven primarily by the tendency to optimize hydrogen bonding [60]. Consistent with the geometry of the S P C / E water molecule, the hydrogen atoms also show a tendency to layer, though it is less pronounced than the oxygens and differs between the K + and C P bearing pores. In particular, there is typically a local maximum in the hydrogen probability density in the K + pores corresponding to hydrogen atoms directed in a more radially outward direction, a tendency also evident in the axial projections in Figs. 7.7 though 7.13. As this radial flaying can be eliminated by reversing the sign of the K + while holding all other parameters constant, the anomaly is likely due to the repulsion between the positive hydrogen atoms and the positive K + ion. Moreover, Chapter 7. W a t e r A d s o r p t i o n i n Ion b e a r i n g Nanopores 121 Figure 7.9: Filling in the 0.54nm pore with a K ion. (a) p = — 66kJ/mol (n = 7) (b) Pw = —64kJ/mol (nw = 9)(c) the low density state at pw = — 60.2kJ/mol (n — 17) (d) the high density state at pw = —60.2kJ/mol (n = 26) (e) Pw = —58.9kJ/mol (n = 31) (f) p = —56kJ/mol (n = 37). + w w w w w an w w examination of the longitudinal projections in Figs. 7.7 through 7.13 indicates that the flaying is more pronounced in the vicinity of the K + ion. In the 0.64 and 0.74nm pores after filling, there is a tendency for the adsorbed water molecules to reside in either a layer coating the inner pore surface or in a shaft of water centred on the axis of the pore. In general, the water molecules in the surface layer participate in more hydrogen bonds than do those in the central region of the pore with, for example, approximately 3 hydrogen bonds per water molecule in the surface layer of the 0.74nm K + pore compared to 2.3 hydrogen bonds per water in the central region of the pore (pw = 56.0kJ/mol). Returning to Fig. 7.3, note that among the cations there is a trend where C a to produce higher water densities than N a , and N a + + 2 + tends tends to produce higher water densities than K . Figure 7.14 shows that this trend is perfectly consistent with the + ion-water interaction energies; with the cations with the stronger ion-water interactions yielding the higher densities in the pore. As the favourable ion-water interaction acts to lower the free energy penalty (approximately equal to —pw) for adsorbing a water Chapter 7. W a t e r A d s o r p t i o n i n Ion b e a r i n g (d) (e) Nanopores 122 (f) Figure 7.10: Filling in the 0.64nm pore with a C P ion. (a) p = -66kJ/mol (n = 6) (b) P w = —64kJ/mol ( n w — 7) (c) the low density state at p w = —60.6kJ/mol (n = 19) (d) the high density state at p = —60.6kJ/mol ( n w = 50) (e) P w = —59.8kJ/mol ( n w = 53) (f) p = —56kJ/mol (n = 63). w w w w w w Chapter 7. Water Adsorption (d) in Ion bearing (e) Nanopores 123 (f) Figure 7.11: Filling in the 0.64nm pore with a K ion. (a) p = -66kJ/mol (n = 6) (b) p = —61.1kJ/mol (n = 14)(c) the low density state at p — -59.9kJ/mol (n = 22) (d) the high density state at p = -59.9kJ/mol (n = 49) (e) p = -58.9kJ/mol (n = 52) (f) p = -56kJ/mol (n = 60). + w w w w w w w w w w w w Chapter 7. W a t e r A d s o r p t i o n i n Ion b e a r i n g Nanopores 124 Figure 7.12: Filling in the 0.74nm pore with a CI ion. (a) pw = —66kJ/mol (nw = 5) (b) p = -64kJ/mol (n = 6) (c) p = -60.3kJ/mol (n = 20) (d) immediately prior to filling pw = —59.9kJ/mol (nw = 25) (e) immediately after filling p = -59.8kJ/mol (n = 83) (f) pw — -56kJ/rnol (n = 93). Parts (a) and (b) show the tendency for water to layer in the wall potential minimum even when only the ion's first solvation shell is present. w w w w w w w Chapter 7. Water Adsorption (d) in Ion bearing (e) Nanopores 125 (f) Figure 7.13: Filling in the 0.74nm tube with a K ion. (a) p = -66kJ/mol (n = 6) (b) p = -64kJ/mol (n = 7)(c) immediately prior to rilling p = -59.7kJ/mol (n = 24) (d) immediately after filling p = -59.4kJ/mol (n = 82) (e) p = -57.5kJ/mol (n = 88) (f) p = -56kJ/mol (n = 92). In (f) a second inner layer of water is beginning to form. + w w w w w w w w w w w w Chapter 7. W a t e r A d s o r p t i o n i n Ion b e a r i n g 10 20 30 40 50 60 15 30 45 Nanopores 60 75 126 90 Figure 7.14: The total water-ion interaction energy (closed points), and the total particlewall interaction energy (open points) for the monovalent ions in (a) 0.44nm, (b) 0.54nm, (c) 0.64nm, and (d) 0.74nm pores. gure 7.15: The total water-water interaction energy (closed points), and the total configurational energy (open points) for the monovalent ions in pores with radius (a) 0.44nm (b) 0.54nm (c) 0.64nm (d) 0.74nm. Chapter 7. Water Adsorption in Ion bearing Nanopores 127 Figure 7.16: Radial probability density profiles for C P pores with radii (a) 0.44nm (b) 0.54nm (c) 0.64nm (d) 0.74nm. In each case, note that the tendency for water to layer the inner surface of the pore is evident even at the lowest chemical potentials. Chapter 7. Water Adsorption in Ion bearing Nanopores 128 Figure 7.17: Radial probability density profiles for K pores with radii (a) 0.44nm (b) 0.54nm (c) 0.64nm (d) 0.74nm. Note that the filled state hydrogen distributions exhibit increased density at large radii compared to the profiles of Fig. 7.16. + Chapter 7. Water Adsorption in Ion bearing Nanopores 129 0.44nm and 0.54nm. In both cases pw = —59.0kJ/mol. molecule, this trend was expected on thermodynamic grounds. Considering the anions, the F ~ ion with the stronger ion-water interactions draws more water into the pore than the C P ion. Comparing the cations and anions however, yields an interesting anomaly with the C P , which typically has the weakest ion-water interactions of all the ionic species except in the 0.74nm pore. Yet, the C P ion produces higher densities than N a and K + + in the pores with radii greater than 0.44nm. To investigate the differences in the filled state densities between the cation and the C P bearing pores, we focus on the behaviour of the C P and K + laden pores with radii of R = 0.44nm and R = 0.54nm. Figure 7.18 shows the regions of the radial distribution functions corresponding to the first solvation shells of the K + and C P ions in the 0.44nm and 0.54nm pores. Here the ion-water radial distribution function is given by Chapter 7. Water Adsorption in Ion bearing Nanopores r~ 5 (Ii 2 130 - r) x (7.22) The first peak in the chloride-water-hydrogen distribution occurs at approximately 0.21nm (pw = -59 kJ/mol) compared to approximately 0.18nm for the pure water peak, while the chloride-water-oxygen peak occurs at approximately 0.31nm. The first potassium-wateroxygen peak occurs at approximately 0.27nm,which is nearly identical to the mean wateroxygen-water-oxygen nearest neighbour distance, while the potassium-water-hydrogen peak occurs at 0.34nm. Figure 7.18 shows that, despite a change in pore radius from 0.44nm to 0.54nm, the typical ion-water distances do not shift significantly. This indicates that the water-ion distances deviate from the typical water-water distances significantly in both pores. The fact that the chloride-hydrogen distance is considerably smaller than the ion-water-oxygen distances is due to the fact that the van der Waals repulsion in the S P C / E water is centred on the oxygens. This allows hydrogen atoms to approach closer to the chloride ion than the oxygen can approach the cations, and gives a more favourable electrostatic interaction than would be expected given the 0.3785nm van der Waals radius of the CP-water Lennard-Jones potential. Despite the fact that K and C P both produce significant deviations in the water + structure in the first solvation shell, the average water-water interaction energies nonetheless show an advantage in the C P carrying pores compared to similar pores with K . + Recall that in the K + bearing pores, the water hydrogens have a larger probability of penetrating into the wall potential than they do in the pores with a C P ion. We can quantify these differences in hydrogen orientation by considering the mean number of hydrogen bonds per water molecule in pores with different ions, when the numbers of adsorbed water molecules are equal. We find that 0.54nm pores with a C P ion, show slightly more hydrogen bonding with an average of 2.7 bonds per water molecule (at Pw = —58.7kJ/mol, with (n ) w — 36.4), compared to 2.4 bonds per water molecule in Chapter the K + pores (at p w 7. Water Adsorption in Ion bearing Nanopores 131 = — 56(kJ/mol), with (nw) = 36.4). As a comparison, each water participates in approximately 3.1 hydrogen bonds in the pure water pores with 37 water molecules. Also note that in the solvation shell of K , water oxygens are pulled towards + the ion, reducing the likelihood of hydrogen bond acceptance from these water molecules. In the C P pore, the images in Fig. 7.8 show that one hydrogen is drawn in toward the C P leaving the other available for hydrogen bonding along with the water oxygen. As a result of these geometric differences, water molecules in the first solvation shell of C P participate in an average of 1.7 hydrogen bonds, while those in the first solvation shell of K + participate in approximately 1.5. These results suggest that in the C P bearing pores, the water structure locally induced by the ion produces a more energetically favourable water geometry than is the case in the K + bearing pores. These differences in hydrogen bonding help explain why the C P bearing pores typically exhibit lower water-water interaction energies than the K + bearing pores, and tend to adsorb a larger number of water molecules in pores with radii greater than 0.54nm. Now consider the 0.44nm pore, where the profiles in Figs. 7.16 and 7.17 show that the C P ion tends to remain in the centre of the pore, while the K + tends to move off axis. If we consider that the C P ion-wall van der Waals radius is 0.38nm, while the radius for K + is 0.33nm, the difference in the mean positions of the two ions can be explained by the size difference between the chloride and the potassium. Given that pure water forms a nearly one dimensional chain in the narrow 0.44nm pore, the fact that chloride is confined to the pore centre explains why there are approximately two water molecules in its first solvation shell. The ability of K + to move off axis while maintaining a favourable interaction with the wall, allows an increase in the number of solvating water molecules. This is shown in the sample configuration displayed in Fig. 7.7(e). Given that the mean difference in total water number between the K + and C P laden 0.44nm pores is approximately one water (at Pw = — 59kJ/mol), the difference in typical solvation number explains the higher water density in the K + 0.44nm pores. For pores larger than 0.54nm, the differences in the filled state densities between the pores with different ions become small compared to the mean water density in the pores. Additionally, the difference in filling chemical potential between the pores decreases. This is the expected behaviour as the ionic concentration Chapter 7. Water Adsorption in Ion bearing drops with increasing R, and in the limit of R 3> o I W Nanopores 132 we expect the isotherms in the ion-laden pores would approach those displayed in Fig. 7.2. 7.3.4 Pressures Having examined the filling isotherms, we now consider the pressures in the pores as the chemical potential is adjusted. Given the cylindrical geometry of the pore, two distinct pressures can be defined. The radial pressure given by (7.23) E E i ( j>i ^ ) + E ^ i ( / 2 - r 0 > i and the axial pressure given by (7.24) In Eqs. 7.23 and 7.24, f\- is the A; component of the force of particle j on particle i, kji is the k component of the vector from j to i, and V = 2ITR L where R is the geometric 2 radius of the pore. In a large system, far from the boundaries we would expect P and L PR to be equal. However, for pores with characteristic radii on the order of a few tenths of a nanometer, we expect the interactions with the bounding wall to cause significant differences in the two pressures. Figure 7.19 displays the radial and axial pressures in K + and C l ~ bearing pores, and shows that after filling the radial pressure typically exceeds the axial pressure by nearly an order of magnitude. The pressures themselves are large, but consistent with pressures computed for pure water pores [60], and pores with an uncharged solute [114]. Note that the first term in Eqs. 7.23 and 7.24 is the ideal gas term which, for the sizes and filled state densities of the pores used here, is on the order of Chapter 4000 • • R= —• R = 3000 — R = *—* R = I 7. W a t e r A d s o r p t i o n i n Ion b e a r i n g 0.44nm 0.54nm 0.64nm 0.74nm Nanopores 133 -1—I—I—I—I—I—I—I- 2000 Figure 7.19: (a) The radial pressures in K (closed symbols) and C l ~ (open symbols) bearing pores (b) the axial pressures. In both cases lines are drawn as a guide to the eye. Errors are approximately 15% of the displayed values + ~ 10 atm. In the case of P , the ideal term is typically larger than the observed pressure, 2 L indicating that the projection of the water-water and water-ion interactions along the pore axis is typically attractive. The large size of the radial pressures can be ascribed to strong repulsive interactions between the wall and the stored water molecules. This is consistent with the observed layering of the adsorbed water displayed in Figs. 7.16 and 7.17 where the water lies within the well of the wall potential, but typically on the steeply repulsive side of the potential minimum. Also note that the axial pressure is negative at values of p corresponding to the filling regions of the isotherms in Fig. 7.3. These w large negative pressures are entirely consistent with the existence of a cluster around the confined ion and indicate that any water molecules located in the ample vacant regions of the pore experience a strong attractive force pulling them toward the cluster. If we now consider that the pressure of saturated water vapour at 298K is approximately 8 x 10~ atm and compare this with the radial pressures of Fig. 7.19(a), it is 3 clear that in order to confine water at the densities displayed in Fig. 7.3 the pore must Chapter 7. Water Adsorption in Ion bearing Nanopores 134 sustain a large pressure gradient. However, nanoscale pores such as single-walled carbon nanotubes are known to have Young's Moduli on the order of ~ 10 GPa [115], and recent 2 studies of the forced insertion of buckminsterfullerene into carbon nanotubes [116] have demonstrated that single-walled carbon nanotubes can withstand pressure gradients up to the GPa range. These studies suggest that pressure gradients on the order of a thousand atmospheres should be sustainable in nanotubes with radii similar to the pores considered here. 7.3.5 Grand Potentials We now address the issue of calculating the grand potentials of the nanoscale pores. When performing a simulation on a bulk system in the grand canonical ensemble, the grand potential is easily calculated from the expression Q = — PV. To see the analogous expressions for the case of water in a pore with smooth walls, consider the differential of the grand potential dfi = -SdT - A-nRLP dR R - 2nR P dL 2 L - n dp (7.25) w where the second and third terms represent the mechanical work done on the pore as a result of radial and axial deformations, respectively. The parameters T and p, w are by definition intensive, while P and P are expected to depend on the radius of the confining R L cylinder. Additionally, note that the wall potential of the cylinder is smooth and thus the cylinder is invariant under arbitrary translations in the z direction. If we hold R constant, and ignore the finite size of the interacting adsorbed particles, this suggests that P and L PR would be independent of L. With these assumptions and under conditions of fixed p,w, T, and R dQ = -2ixR P dL. (7.26) 2 L As Pi is independent of L we can trivially integrate to find that tt(p , T, R, L) = -2nR LP 2 w L = -P V. L (7.27) Chapter 7. Water Adsorption in Ion bearing Nanopores 135 An alternative route to this result is to note that, with the assumptions given above, we would expect the total internal energy U to behave like a first order homogeneous function of the extensive variables 5, V, L, and nw Assuming this holds for arbitrary scalings of the extensive variables, then introducing the scale factor n allows us to write nU(S, V,L,n ) w = U(nS,r]V,riL,r]n ) w (7.28) and by Euler's theorem for homogeneous functions [117], U = TS — PRV — (Pi — PR)V + Hw^w- Legendre transforming into the grand canonical ensemble again gives Q — —PiV. However, in the pores we consider in this chapter, the total pore length is only about nine times the van der Waals radius of a water molecule. As a result it is not clear that the arguments for the intensive behaviour of Pi given above are applicable. This is particularly true for the ion-bearing pores where scaling L changes the ionic concentration in the pore. With the loss of intensivity of the axial pressure, trivial integration of Eq. 7.26 cannot be performed. In the limit of large L however where the effects of the finite water size or the single ion in a period would be small, we would again expect Eq. 7.27 to be a valid expression for the grand potential. More formally, let P be the pressure in the pore in the limit of large L, and let e be a length where P (L L = e)-P P < 1. (7.29) Let QQ = —PV and then consider the quantity (7.30) Now let 7 = max (\P (x)\) L on 0 < x < e, then Chapter 7. Water Adsorption in Ion bearing Nanopores 136 L (nm) Figure 7.20: The behaviour of the axial pressure PL under scaling of the pore length L. While the pressure in the pure water pore oscillates as the length is scaled, the variation in Pi over the range of displayed values of L is much smaller than the variation in the K bearing pore (pw = —57.0kJ/mol). + < (7.31) 1 + P Equation 7.31 shows that the trivial integration of Eq. 7.26 will only provide reasonable grand potentials in the limit that L ^> e. In other words, and consistent with our expectation, the trivial integration of the grand potential is only valid in the thermodynamic limit of large L. Now consider Figure 7.20 which shows P L in a K + as a function L in a pure water pore and bearing pore both with R = 0.44nm, and p-w = —57kJ/mol. Note that while the pressure in the pore containing pure water oscillates at small values of L the value of P L averaged over these oscillations is close to the value of P L at large L. Conversely, the ion bearing pore shows a large nearly monotonic decrease in P L The results of Fig. 7.20 suggest that for a K + as L is increased. bearing pore e in Eq. 7.31 is « 3.0nm. This indicates that for ion-bearing pores with L = 1.5nm, Eq. 7.27 would yield a poor approximation to the grand potential. In the pure water pore, it would be appropriate Chapter 7. Water Adsorption in Ion bearing Nanopores 137 to place e in Eq. 7.31 between 1.5nm and 2.0nm, and thus while Eq. 7.27 would not be expected to provide an accurate estimate of the grand potential, it would yield better results for the pure water case than for the ion-bearing pore. Also note that if the average over the oscillations observed in Pi for pure water produces a value of the pressure close to the limiting value of Pi, Eq. 7.27 will perform better than would be expected based on the error estimate in Eq. 7.31. As an alternative to Eq. 7.27, one possible means of computing accurate grand potentials is to integrate the curves in Fig. 7.20. However, the large errors in the pressures make this route to the grand potential untenable. Another means of calculating Q, is to use the adsorption isotherms we have already computed for these pore as follows (7.32) In Eq. 7.32, p w is a very low value of the chemical potential at which (nw) ~ 0. At this low chemical potential, the only particle occupying the pore is the ion and the nw = 0 term of Eq. 7.6 can be readily integrated giving Q(p \v, l w discussion that follows, pff L, T) = — k Tm(Zo). B In the = — 115kJ/mol unless otherwise noted, and all grand poten- tials are computed using the composite Simpson's rule [3] with a step-size of 0.15kJ/mol. Figure 7.21 shows the grand potential of K + and C P bearing pores of varying radii as a function of the reservoir chemical potential. In all cases the grand potential of the K bearing pores was found to be lower than the C P pores indicating that the K + + pores are more thermodynamically stable. The difference is most significant in the smallest pores, and in the larger pores the difference is reduced as the chemical potential is increased and the pore fills. The reduction in the grand potential difference on filling is consistent with the earlier conclusion that in the filled state of the larger, R = 0.64nm and R = 0.74nm pores, the water dominates the thermodynamics of the pore with the ion becoming a comparatively minor effect. Figure 7.22 shows the per particle entropy in each pore derived from the grand potentials, and demonstrates that in pores with R > 0.54nm the per particle entropy is Chapter 7. W a t e r A d s o r p t i o n in Ion b e a r i n g Nanopores 138 -100 -125 -150 -175 -200 E -63 -62 -61 -60 -59 -58 """-63 -62 -61 -60 -59 -58 G -150 -200 -250 -300 ~ -63 3 5 -62 -61 -60 -59 -58 H w -63 -62 -61 -60 -59 -58 (kJ/mol) Figure 7.21: The grand potential of K (circles) and C P (squares) bearing pores with (a) R = 0.44nm, (b) R = 0.54nm, (c) R = 0.64nm, and (d) R = 0.74nm. Energies are in (kJ/mol) of pore. + higher in K + bearing pores than in C P bearing pores given the same number of adsorbed water molecules. This is consistent with the higher mean number of hydrogen bonds per water molecule observed in C P bearing pores, and the more favourable water-water interaction energies displayed by the pore with a C P ion. In the R = 0.44nm pore, a crossover is observed in the per particle entropy with the C P pore showing the larger per particle entropies at low chemical potentials. This could be predicted from the energetic trends discussed previously where the water-K interaction is stronger than the water-CP + interaction. Therefore at low chemical potentials, where there are relatively few water molecules, the behaviour of the water will be dominated by the ion-water interaction. The stronger interaction between the K + ion and the water would then be expected to lead to lower per particle entropies as the water is more tightly held about the ion. Conversely, as the number of adsorbed water molecules increases and the mean interaction strength between the ions and the adsorbed water molecules decreases, the stronger water-water interactions that occur in the C P laden pore would be expected to yield the lower per particle entropy. Note that while not shown, in the pores with R > 0.54nm, a similar crossover in the per particle entropies occurs at low chemical potentials. Chapter 7. Water Adsorption in Ion bearing Nanopores 139 0.08 0.076 0.072 * — s | I 6 8 10 12 14 10 15 20 25 30 35 0.09 n 0.084 0.078 0.072 0.066 °'°^0 20 30 40 50 6 0 2 0 40 60 8fT~ Figure 7.22: The entropy in K (circles) and C P (squares) bearing pores with (a) R = 0.44nm, (b) R = 0.54nm, (c) R = 0.64nm, and (d) R = 0.74nm. Entropies are in (kJ/mol • K) of particles. + 7.3.6 Boundary Effects We now address the topic raised earlier about the influence of the periodic boundary conditions on our results. Beginning with the issue of the polarization interface induced by the periodic boundary conditions, we note that the interfacial energy should scale roughly as ~ R . 2 We reach this conclusion by assuming that the interactions between water molecules are short-ranged which would imply that the interfacial energy should scale with the cross sectional area of the pore. The assumption of a short-ranged waterwater interaction can be justified by noting that the water molecules can be approximated as dipoles with moments of approximately 7.8 x 10~ C • m. Setting r to be the distance 30 between two water molecules, the water-water interaction energy would then scale as r ~ if 3 the relative orientation of the molecules were fixed, and as r ~ if there is thermal averaging 6 over the orientations. If the interfacial energy goes as ~ R , then if R is fixed and L is 2 scaled, the interfacial energy will decrease as a fraction of the total energy in the system at least as fast as ~ L~ , assuming the total energy in the filled state is proportional to l the total volume. Based on these scaling arguments, we expect that as L is increased holding R fixed Chapter 7. W a t e r A d s o r p t i o n i n Ion b e a r i n g Nanopores 140 •— L= 1.50 n m 80 *—L = 3.10nm 60 40 20 : 0 . i -64 -66 ( a ) . . . > . . , -62 i . . . . i . . , . ) , . . , -60 -58 -56 i 1 1 1 1 (d) 1 r-54 1 Hydrogen A • 0 0.1 0.2 0 0.1 0.2 o0 1111111 0.2 0.1 0.3 0.4 r (nm) Figure 7.23: (a) The adsorption isotherms for K bearing pores with radius R = 0.54nm with different lengths. The (b) ion, (c) water oxygen, and (d) water hydrogen radial probability profiles. The points are for the L = 1.5nm pore while the lines are for the L = 3.1nm pore (p = —58.2 kJ/mol). + w whatever influence the interface induced by the periodic boundaries has on the behaviour of the system will decrease as roughly ~ L _ 1 , or in other words at roughly the same rate as the ion density. Figure 7.23(a) shows isotherms for two R = 0.54nm pores, one with L — 1.5nm, the other with L = 3.1nm, and each containing one K + ion. The plot shows that increasing L has the effect of sharpening the filling transition and increasing the filling value of p\y. If we consider the values of (nw) immediately after filling, we find that (nw(L = 1.5)) ~ 33, (nw(L = 3.1)) ~ 69, which is nearly perfect linear scaling with L. In going from L = 1.5nm to L — 3.1nm, the total energy changes by a factor of approximately 2, as does the wall energy. The magnitude of the water-solute energy only increases by approximately 10%. However, the magnitude of the water-water interaction energy increases by a factor of approximately 2.4. This super-linear change in the waterwater interaction energy can be ascribed both to a reduction in the boundary effect and the fact that there is less disruption of the water structure in regions of the pore between the ion and the boundaries of the simulation cell. Figure 7.23(b) through (d) displays the radial probability density for the various species in the two pores. There clearly is Chapter 7. W a t e r A d s o r p t i o n i n Ion b e a r i n g K Nanopores 141 Periodic • • K Capped • * K° Capped + (b) Oxygen 0.1 0.2 r (nm) 0 0.15 0.3 0.45 Figure 7.24: The (a) solute, (b) water oxygen, and (c) water hydrogen radial probability profiles for the different species in capped and periodic pores. K° is a neutral solute with the same Lennard-Jones parameters as K . In all cases the pore radius was R = 0.28nm and the pore was filled with nw = 34 water molecules. + very little change in going from a L = 1.5nm pore to a L = 3.1nm pore. The reduction in the secondary peak in the hydrogen distribution displayed in Fig. 7.23(d), is expected as the distribution represents an average over the length of the pore. As the secondary hydrogen peak is due to the repulsion between the positive K + ion and the hydrogens, and is strongest close to the ion, lengthening the pore would be expected to reduce the significance of this effect in the average. The results of Fig. 7.23 indicate that when L — 1.5nm, the simulation results for the water structure and ion solvation in the filled pores are not dominated by effects introduced by the boundary conditions. Returning to the discussion in the introduction to this chapter, this means that it is not unreasonable to view the results in this chapter as being an approximation to, for example, a finite 3.0nm channel with open ends. In this interpretation, the ion is confined to the centre of the pore in the direction parallel to the pore axis, while being free to move in the plane of the cylinder. To further explore the influence of boundary conditions on the water structure and ion solvation, we can consider the effects of capping the ends of the pores. Here we Chapter 7. Water Adsorption in Ion bearing Nanopores 142 consider two pores, where the first uses wall potential of Eq. 7.1 with periodic boundary conditions in the z direction and the minimum image convention. The second pore uses Eq. 7.1 to model the cylindrical wall, supplemented by completely hydrophobic walls that are perpendicular to the pore axis and are placed at z = ±L. These end caps are modelled using the exponential repulsion given in Eq. 7.1, yielding the following expression for the end-cap potential 0 = 2^ e-' cosh(/^). (7.33) L c a p 0 In Eq. 7.33, 4>Q and I have the same values as in Eq. 7.1. Figure 7.24 shows the radial probability distributions for a periodic pore with L = 1.5nm containing a K ion, a capped + pore with L = 1.52nm containing a K + ion, as well as a capped pore with L — 1.52nm containing an uncharged solute particle with the same Lennard-Jones parameters as the K . In all cases the geometric radii of the pores were set to 0.28nm, yielding physical radii + of approximately 0.26nm. Additionally, in all cases the pores were filled with nw = 34 water molecules. The capped pores are slightly longer than the periodic pores to account for the finite range potential of the capped ends, thus ensuring the volume of the pores is similar in the capped and periodic cases. The water oxygen, water hydrogen and K + profiles displayed in Fig. 7.24 show that the switch from periodic to capped ends has only a minor effect on the average water structure and ion radial position in the pores. A comparison of the water-water and water-ion interaction energies between the capped and uncapped pores shows differences of less than 10% in both quantities. The small differences in energy between these two quantities indicate that the solvation of the ion in the capped pores is essentially unchanged from the periodic pore. A n analysis of the axial probability distribution for the ion in the capped pore shows that the ion is rarely less than approximately 0.3nm from one of the pore ends. In pores using Eq. 7.2 with a radius of 0.54nm, the typical distance between the K + ion and a water oxygen in the ion's first solvation shell is approximately 0.27nm. Assuming the radius of the first solvation shell in the capped pores is also approximately 0.27nm, we can conclude that the ion rarely moves close enough to the end caps to strip away the first solvation shell, a behaviour that would be expected with a highly hydrophilic solute like K . Conversely, an analysis + Chapter 7. Water Adsorption in Ion bearing Nanopores 143 of the behaviour of an uncharged solute in the capped pores shows that, as expected, the neutral solute is typically excluded from the water by being forced to one of the ends of the pores. This is the cause of the significant differences in the solute profiles shown in Fig. 7.24(a). These results indicate that by far the largest factor affecting the filled state structures in the pores is the cylindrical boundary. Even in a relatively short 3.0nm pore, changing the nature of the boundary conditions, as well as the solute does not change the qualitative character of the water structure in the pores. 7.4 Summary Through the use of grand canonical Monte Carlo simulations we have examined the adsorption of water into carbon-nanotube-like pores 3.0nm in length containing single ions. At reservoir chemical potentials below approximately —62.0kJ/mol, water was present in the pores in the form of a cluster forming a solvating shell around the ion. This water-ion cluster tends to lie near the minimum of the wall potential. In pores with radii of 0.44nm and 0.54nm, increasing the chemical potential resulted in gradual growth of the water cluster around the ion until the pore was effectively filled with water. In larger 0.64nm and 0.74nm pores, the filling process was more abrupt with gradual growth of the ion cluster to the point where the cluster covered the cross section of the pore, followed by rapid filling of the entire pore on an increase in the chemical potential. Regardless of size and ionic species, the pores filled at chemical potentials lower than the bulk chemical potential of water at 298K. If we consider the narrow pores to be analogous to biological ion channels, these results suggest the presence of a single ion in the channel will likely cause the channel to flood, assuming the channel has time to come to equilibrium before the ion traverses the channel. Note that as the pore radius increased, and hence the ion density decreased, the filling chemical potential showed the expected trend and increased as well. In the filled state of the smallest 0.44nm pore, the water molecules tend to form a line with a highly ordered dipole alignment in the directions favoured by the ion. In the 0.54nm pore, the water forms a hydrogen-bonded structure with a cross section that can be approximated by a square with a side length of approximately 0.27nm. In the larger Chapter 7. Water Adsorption in Ion bearing Nanopores 144 pores, the water forms a layer that coats the inner surface of the pore, as well as a line of water down the centre of the pore in the 0.64nm case. At high chemical potentials, the 0.74nm pore exhibits the growth of a second inner water layer. In general, water in filled ion-free pores and filled ion-laden pores of the same radius forms very similar structures, indicating that cylindrical confinement is the primary source of the observed structure. The presence of the ion does have the effect of ordering the water dipole moments, particularly in the 0.44nm pore, and causes a local disruption of the water structure in the region about the ion. This is particularly evident in the pores with cations where the water molecules near the ion tend to orient with their hydrogens radially outward, away from the ion. This results in a noticeable increase in the radial probability distribution for water hydrogens at large radii in cation bearing pores. In all pore sizes and with all ionic species, the density of water in the filled state of the pores exceeded the density of liquid water at 298K. While the filling behaviour of the nanopores was qualitatively the same regardless of ionic species, comparison of C P bearing pores to K + bearing pores showed more favourable hydrogen bonding in the C P pores. In larger pores, where the water-water interactions dominate, the more stable water structure associated with the negative chloride ion results in a higher water density in the filled state. Conversely, comparing pores with R = 0.44nm, K + bearing pores showed a larger water density. This was shown to be the result of the smaller size of the K + ion. Grand potentials were computed for pores with C P and those with K , and despite the higher water densities in the C P pores, the + grand potential was found to be lower in the K + case. Per particle entropies were also found to be lower in the C P pore consistent with the finding of a more cohesive water structure in pores with this ion. In all cases the radial pressures in the pores in the filled state were found to be on the order of 10 atm. 3 At chemical potentials corresponding to the filled states in the ion-bearing pores, saturated water vapour is well approximated by an ideal gas with pressures well below one atmosphere. While the resulting pressure gradients across the pore are large, structural studies of nanoscale cylinders indicate that single-walled carbon nanotubes with radii comparable to those investigated in this chapter, should be able Chapter 7. Water Adsorption in Ion bearing Nanopores 145 to sustain gradients of this magnitude. This suggests that carbon nanotubes are viable structures for nanoscale fluid containment devices. 146 Chapter 8 Conclusion and Suggestions for Future Work This thesis has been concerned with the development of models for, and an examination of, transport and structure formation in nanoscale channels. The first set of results presented were based on the Totally Asymmetric Simple Exclusion Process (TASEP), a stochastic model that is suitable for describing the unidirectional transport of discrete objects along one dimensional pathways, when the transported objects move in discrete steps. In Chapter 3, an extension to this model was introduced where the size of the transported particles exceeded the size of the steps the particles make by a factor of d. Monte Carlo simulations of this modified T A S E P showed that the T A S E P with extended particles retained the three-phase character of the standard T A S E P . By explicitly enumerating the possible particle occupancy states in the interior of the lattice, an accurate mean field estimate for the current and the bulk density in the maximal current phase was introduced. Similarly, by explicitly enumerating the possible occupancy states of the first and last d sites of the lattice, accurate current and density estimates for the entry limited phases were also developed. This extension of the T A S E P allows it to be applied directly to a variety of biophysical problems such as ribosome and molecular motor motion [34]. Although in Chapter 3 we considered only total asymmetry, in future work analogous expressions could be computed for partially asymmetric models by allowing backward particle steps and extraction and injection from both ends of the T A S E P lattice. Another possible variation is to consider TASEPs where particles of different sizes are allowed to occupy the lattice simultaneously. For example, in the maximal current phase the partition function analogous to Eq. 3.11 would be Z({ },L) ni (Eini)!/iV-Ei"i(rfi-l) ) (8.1) Chapter where 8. Conclusions and Suggestions for Future Work 147 is the number of particles of size i in an L-site lattice segment. This expression yields the following approximation for total particle number current (8.2) where Oj is the bulk density of particles of size j. The current J could then be maximized with respect to the o given specific constraints such as Oi/o x where a m i+1 = ai/a i i+ in the case < j3 for all m and n. This preliminary analysis suggests that there are no n regimes in which J would depend only on the interior hopping rates Pi. The analogous boundary limited states are complex given the possible dependence on the injection and extraction rates for all the different particle species, and would require extensive Monte Carlo simulations to investigate. Continuing with modifications to the TASEP, Chapter 4 considered the case where there are defect lattice sites arranged periodically with period T. Any particle that occupies a defect lattice site moves to the right at a reduced rate p 2 Monte Carlo simulations were again used to explore the phase diagram of this modified TASEP, which retains the three phase character of the standard TASEP. Analytic expressions for the bulk current and density in the three phases of the T A S E P were developed for the T = 2 case, and the finite segment mean field (FSMF) method was developed to provide numerical approximations to the current and bulk densities for T > 2. In the T A S E P phase diagram, the modifications associated with the introduction of the second hopping rate were an increase in the size of the maximal current phase as p was reduced, and a loss of symmetry 2 of the phase diagram about the a = f3 line for certain values of of the lattice size N, the defect period T , and the location of the first defect 5. In future work, an exploration of the correlations within a single period and between periods of the periodic T A S E P could be performed, both for intrinsic interest in the behaviour of this variant of the T A S E P and in hopes of further refining the finite-segment mean field (FSMF) method. Additionally, note that in the work presented here the spatial variation in the movement rates was sharp; the lattice was perfectly uniform except for a finite number of slow defect sites. A n interesting extension would be to consider the case where the spatial variation in the hopping rates was smooth, with a gradual change in Chapter 8. Conclusions and Suggestions for Future Work 148 the hopping rate extended over a large number of lattice sites. In this case Eq. 2.4 could be used to derive a mean field, continuum, equation [118, 119] for the densities along the entire length of the T A S E P lattice. The solution to this equation would provide mean field estimates for the current, and the entire particle density profile. Additionally, the solution to the continuum equations could be used to study how the quality of the mean field estimates depend on the functional form of the spatially varying hopping rates. In Chapter 5, the T A S E P and the finite segment method developed in Chapter 4 were used to study ribosome movement along messenger RNAs that include slowly translated codons. Not unexpectedly, groups of slow codons were found to produce more of a reduction in the ribosome current than isolated codons. However, the additional reduction in the ribosome current gained by increasing the cluster size of the slow codons beyond approximately five codons was found to be minimal. When two successive slow codons were spaced more than eight normal codons apart, the effect of the second codon was found to be minimal and the ribosome current was approximately equal to the current when only a single slow codon was present. Both the requirement for close spacing for synergistic effects between multiple codons, and the approximately five codon limit on the slow codon cluster size, match well with the experimental observation of the grouping of slow codons in the mRNAs of E. coli bacteria. We can speculate that the existence of clusters of slow codons in real mRNA may have the role of providing pause points for partial, perhaps chaperone assisted, folding of the growing polypeptide. If this is the case, the experimental data on slow codon frequency could be expected to show significant separations between slow codon clusters in order to ensure that the clusters provided only local pause points, while producing only a minimal overall reduction in the ribosome current. The dehybridization and translocation of D N A through nanopores is another biologically related transport process, and was investigated in Chapter 6. Following recent experimental work, we considered the voltage driven translocation of D N A hairpins through membrane embedded a-hemolysin channels, using a two dimensional stochastic model that coupled the dehybridization and sliding of the hairpin. In the limit of low voltage, the D N A hairpin was found to escape the a-hemolysin channel by sliding in a direction Chapter 8. Conclusions and Suggestions for Future opposite to the one favoured by the applied voltage. Work 149 As the voltage was increased, a distinct peak in the mean escape time was predicted. As yet, experiments have not been performed to verify the existence of this predicted peak in the escape time. The peak occurred at a voltage where the characteristic time scale for reverse escape was comparable to the characteristic time scale for dehybridization of the hairpin. In future work, many improvements and modifications to the model could be made. For example, the estimates for the entropic penalty associated with the dangling double stranded D N A segment in the cis chamber could be refined. Additionally, an estimate for the value of UR, the rate of rebinding of a completely dehybridized hairpin, could be made by considering the behaviour of a polymer confined in the vestibule region. If dehybridized D N A bases can only rehybridize when they pass within a certain minimum distance of each other, an investigation of the polymer dynamics of the dehybridized D N A hairpin could yield an approximate time scale for rehybridization of any two complimentary bases. Finally, beyond all the applications to D N A translocation, the model in Chapter 6 suggests an interesting problem in the theory of random walks. The two dimensional model presented in Chapter 6 can be reduced to the problem of two random walkers on a finite one-dimensional lattice with absorbing boundaries at each end. The random walkers move independently of each other, with spatially varying stepping rates, subject to the condition that the walkers can never cross. The question is to calculate the mean time required for one of the walkers to reach one of the absorbing boundaries. The analogous problem of a random walk to an absorbing state on a lattice with quenched disorder has been solved exactly [120] for a single walker, and the two walker problem is a natural and interesting extension. Turning to the subject of structure formation in nanoscale channels, in Chapter 7 grand canonical Monte Carlo simulations were used to investigate the behaviour of water adsorbed in ion-bearing nanoscale pores. The simulations showed that at low water reservoir chemical potentials, a small number of water molecules were adsorbed in a cluster around the single ion in the pores. As the chemical potential was increased past approximately — 59kJ/mol, all pores were considered to be filled with water. In the smaller 0.44nm and 0.54nm pores, the filling occurred through a gradual growth of the Chapter 8. Conclusions and Suggestions for Future Work 150 water cluster found around the ion. Conversely in the 0.64nm and 0.74nm pores, the water clusters were observed to grow until the cluster covered the pore cross section, a small increase in chemical potential would then completely fill the pore. The structures formed by the water in the filled states of the pore depended strongly on the pore radius. In the 0.44nm pore, water was found to form a nearly one-dimensional hydrogen-bonded line with the water dipoles oriented in the direction favoured by the ion. In the 0.54nm pore, a hydrogen bonded structure was formed that had an approximately square cross section with sides of approximately 0.27nm. In the 0.64 and 0.74nm pores, the water formed a layer one water molecule thick, coating the inside of the cylindrical wall. Additionally, in the 0.64nm pore a shaft of water arranged in nearly a line was found in the centre of the pore, while in the 0.74nm pore a second cylindrical layer of water partially formed. In all cases the density of water in the ion-bearing pores after filling was found to exceed the density of bulk liquid water at 298K and latm pressure. However, in the filled state the water density in ion-bearing and ion-free pores of the same radii were comparable. Moreover, the structures formed by water adsorbed in ion-free and ion-bearing pores are qualitatively similar except in the region near the ion. These observations suggest that the major effect of the ion is to reduce the filling chemical potential of the pore, without globally disrupting the filled-state water structures. The pressures in the water-filled pores were also investigated and the radial pressures were found to be on the order of 10 3 atm. This suggests that under ambient conditions, a filled nanoscale pore would need to sustain a significant pressure gradient. While large, these gradients have been shown to be in the range of tolerance for nanostructures such as single-walled carbon nanotubes. In terms of future work, changing the boundary conditions by exposing the ends of the pore to reservoirs of bulk water would yield a more realistic simulation of systems such as biological ion channels. Another interesting extension would be to allow the ions themselves to be transfered to and from the particle bath. This would more realistically model the behaviour of a nanopore or ion channel immersed in, for example, a salt solution. Early investigations have shown that while ions will readily enter the pore, it is very difficult to remove them due to the large energy barrier associated with breaking anioncation pairs, or pulling an ion out of a closely-packed water cluster. Without ion extraction Chapter 8. Conclusions and Suggestions for Future Work 151 transitions, the Monte Carlo simulation does a poor job of sampling the space of possible configurations, and the Monte Carlo results are suspect. One way around this difficulty would be to devise a means of inserting and extracting ion-water clusters rather than single ions. In addition to changing the boundary condition at the ends of the pore, another possible extension to this research would be to change the wall potential by decorating it with charged groups. In biological ion channels, charged carboxyl and amino groups are typically found near the channel extremities, and are known to play a critical role in the observed selectivity of the channels for specific species of ions. However, it would also be interesting to observe the influence of these groups on the structure of absorbed water. In particular, the nucleation of water clusters near the charged groups would create solvation zones where the water-ion interaction would be favourable. The gaps between these zones, being empty, would represent a potential barrier the ion would have to traverse in order to move from one side of the channel to the other. As the strength of these potential barriers would affect the rate at which ions traversed the channel, an understanding of how the strength and arrangement of the charged groups affects the adsorbed water density and structure would aid in the development of more coarse-grained stochastic models of ion transport. 152 Bibliography [1] R. B. Lehoucq, D. C. Sorensen, and C. Yang, A R P A C K Users' Guide: t i o n of L a r g e - S c a l e E i g e n v a l u e Problems with Implicityly Restarted A r n o l d i (SIAM, Philadelphia, Pennsylvania, 1998). SoluMethods [2] M A T L A B , The MathWorks Inc., 3 Apple Hill Drive Natick, M A 01760-2098 (2005). [3] W. H. Press, S. A . Teukolsky, W. T . Vetterling, and B. P. Flannery, N u m e r i c a l Recipies i n C (Cambridge University Press, Cambridge, United Kingdom, 1995), 2nd ed. [4] M . Zuker, Nucleic Acids Research 31, 3406 (2003). [5] J. A . Spudich, Nature 372, 515 (1994). [6] M . Schliwa and G . Woehlke, Nature 422, 759 (2003). [7] B. Alberts, D. Bray, J. Lewis, M . Raff, K . Roberts, and J. D. Watson, M o l e c u l a r Biology of the C e l l (Garland Publishing, New York, New York, 1994), 3rd ed. [8] J. Solomovici, T . Lesnik, and C. Reiss, Journal of Theoretical Biology 1 8 5 , 511 (1996). [9] T . J. Minehardt, N. Marzari, R. Cooke, E . Pate, P. A . Kollman, and R. Car, Biophysical Journal 8 2 , 660 (2002). [10] J. D. Lawson, E . Pate, I. Rayment, and R. G . Yount, Biophysical Journal 8 6 , 3794 (2004). [11] C. T . MacDonald, J. H . Gibbs, and A . C. Pipkin, Biopolymers 6, 1 (1968). [12] C. T . MacDonald and J. H. Gibbs, Biopolymers 7, 707 (1969). [13] B. Derrida, E . Domany, and D. Mukamel, Journal of Statistical Physics 6 9 , 667 (1992) . [14] S. A. Janowsky and J. L . Lebowitz, Physical Review A 4 5 , 618 (1992). [15] B. Derrida, M . R. Evans, V . Hakim, and V. Pasquier, Journal of Physics A 26, 1493 (1993) . [16] G . M . Schiitz and E . Domany, Journal of Statistical Physics 7 2 , 277 (1993). [17] A . B. Kolomeisky, Journal of Physics A 3 1 , 1153 (1998). Bibliography 153 [18] T . Nagatani, Journal of Physics A 2 6 , 6625 (1993). [19] S. Sandow, Physical Review E 5 0 , 2660 (1994). [20] M . Schreckenberg, A . Schadschneider, K. Nagle, and N. Ito, Physical Review E 5 1 , 2939 (1995). [21] F . H. Essler and V . Rittenberg, Journal of Physics A 2 9 , 3375 (1996). [22] B. Derrida and M. R. Evans, in N o n e q u i l i b r i u m S t a t i s t i c a l M e c h a n i c s i n One D i m e n sion, edited by V . Privman (Cambridge University Press, Cambridge U K , 1997). [23] G. Tripathy and M . Barma, Physical Review Letters 7 8 , 3039 (1997). [24] T . Sasamoto and M . Wadati, Journal of Physics A 3 1 , 6057 (1998). [25] G. Tripathy and M . Barma, Physical Review E 5 8 , 1911 (1998). [26] B. Derrida, Physics Reports 301, 65 (1998). [27] A . B. Kolomeisky, G . M . Schiitz, and E . B. Kolomeisky, Journal of Physics A 3 1 , 6911 (1998). [28] M . Bengrine, A . Benyoussef, H. Ez-Zahraouy, J . Krug, M . loulidi, and F . Mhirech, Journal of Physics A 3 2 , 2527 (1999). [29] V . Karimipour, Physical Review E 5 9 , 205 (1999). [30] V . Popkov and G . M . Schiitz, Europhysics Letters 4 8 , 257 (1999). [31] K. M . Kolwankar and A . Punnoose, Physical Review E 6 1 , 2453 (2000). [32] Z. Nagy, C. Appert, and L. Santen, Journal of Statistical Physics 109, 623 (2002). [33] T . Chou, Biophysical Journal 8 5 , 755 (2003). [34] L. B. Shaw, R. K . P. Zia, and K. H. Lee, Physical Review E 6 8 , 021910 (2003). [35] A . John, A. Schadschneider, D. Chowdhury, and K. Nishinari, Journal of Theoretical Biology 231, 279 (2004). [36] L. Stryer, B i o c h e m i s t r y (W. H. Freeman and Co., New York, New York, 1995), 4th ed. [37] S. E . Henrickson, M . Misakian, B. Robertson, and J . J. Kasianowicz, Physical Review Letters 8 5 , 3057 (2000). [38] A. Meller, L . Nivon, and D. Branton, Physical Review Letters 8 6 , 3435 (2001). [39] S. Howorka, S. Cheley, and H. Bayley, Nature Biotechnology 1 9 , 636 (2001). [40] W. Vercoutere, S. Winters-Hill, H. Olsen, D. Deamer, D. Haussler, and M . Akeson, Nature Biotechnology 1 9 , 248 (2001). Bibliography 154 [41] M . Bates, M . Burns, and A . Meller, Biophysical Journal 84, 2366 (2003). [42] W. Sung and P. J. Park, Physical Review Letters 77, 783 (1996). [43] D. K . Lubensky and D. R. Nelson, Biophysical Journal 77, 1824 (1999). [44] T . Ambjornsson, S. P. Apell, Z. Konkoli, E . Di Marzio, and J. J . Kasianowicz, Journal of Chemical Physics 117, 4063 (2002). [45] C. Y . Kong and M . Muthukumar, Electrophoresis 23, 2697 (2002). [46] E . Slonkina and A . B. Kolomeisky, Journal of Chemical Physics 118, 7112 (2003). [47] O. Flomenbom and J . Klafter, Physical Review E 68, 041910 (2003). [48] S. Cocco, R. Monasson, and J. F . Marko, Proceedings of the National Academy of Sciences of the United States of America 98, 8608 (2001). [49] D. K . Lubensky and D. R. Nelson, Physical Review E 65, 031917 (2002). [50] S. M . Bhattacharjee and D. Marenduzzo, Journal of Physics A 35, 349 (2002). [51] U. Gerland, R. Bundschuh, and T . Hwa, Physical Biology 1, 19 (2004). [52] J. Mathe, H. Visram, V . Viasnoff, Y . Rabin, and A. Meller, Biophysical Journal 87, 3205 (2004). [53] J . Nakane, M . Wiggin, and A . Marziali, Biophysical Journal 87, 615 (2004). [54] W. Nonner, D. Gillespie, D. Henderson, and B. Eisenberg, Journal of Physical Chemistry B 105, 6427 (2001). [55] T . W. Allen, O. S. Andersen, and B. Roux, Proceedings of the National Academy of Sciences of the United States of America 101, 117 (2004). [56] S. Y . Noskov, S. Berneche, and B. Roux, Nature 431, 830 (2004). [57] M . Carillo-Tripp, H . Saint-Martin, and I. Ortega-Blake, Physical Review Letters 93, 168104 (2004). [58] Z. Yang, T . A . V . der Straaten, U . Ravaioli, and Y . Liu, Journal of Computational Electronics 4, 167 (2005). [59] C. Peter and G. Hummer, Biophysical Journal 89, 2222 (2005). [60] R. M . Lynden-Bell and J. C. Rasaiah, Journal of Chemical Physics 105, 9266 (1996). [61] A . Waghe, J. C. Rasaiah, and G. Hummer, Journal of Chemical Physics 117, 10789 (2002). [62] J. Wang, Y . Zhu, J . Zhou, and X . H. Lu, Physical Chemistry Chemical Physics 6, 829 (2004). Bibliography 155 [63] B. D. Huang, Y . Y . Xia, M . W. Zhao, F . L i , X . D. Lui, Y . J. Ji, and C. Song, Journal of Chemical Physics 8, 084708 (2005). [64] B. M . Kim, S. Qian, and H. H. Bau, Nano Letters 5, 873 (2005). [65] D. Chowdhury, L. Santen, and A. Schadschneider, Physics Reports 329, 199 (2000). [66] N. Metropolis and S. Ulam, Journal of the American Statistical Association 44, 335 (1949). [67] A . B. Bortz, M . H . Kalos, and J. L. Lebowitz, Journal of Computational Physics 17, 10 (1975). [68] M . E . J. Newman and G . T . Barkema, M o n t e C a r l o Methods (Oxford University Press, Oxford, United Kingdom, 1999). i nStatistical Physics [69] L. B. Shaw, J . P. Sethna, and K. H. Lee, Physical Review E 7 0 , 021901 (2004). [70] D. Chowdhury and J . S. Wang, Physical Review E 6 5 , 046126 (2002). [71] P. Garrido, J . Marro, and R. Dickman, Annals of Physics 199, 366 (1990). [72] G. Szabo, A. Szolnoki, and L. Bodocs, Physical Review A 4 4 , 6375 (1991). [73] M . Bulmer, Genetics 129, 897 (1991). [74] G . Marais and L . Duret, Journal of Molecular Evolution 5 2 , 275 (2001). [75] J. R. Powell and E . N. Moriyama, Proceedings of the National Academy of Sciences of the United States of America 9 4 , 7784 (1997). [76] C. G. Kurland, F E B S Letters 285, 165 (1991). [77] A . S. Spirin, C e l l - F r e e T r a n s l a t i o n Systems 2003). (Springer-Verlag, New York, New York, [78] M . Robinson, R. Lilley, S. Little, J. S. Emtage, and G . Yamamoto, Nucleic Acids Research 1 2 , 6663 (1984). [79] M . A. Sorensen, C. G . Kurland, and S. Pedersen, Journal of Molecular Biology 365 (1989). [80] D. A . Phoenix and E . Korotkov, F E M S Microbiology Letters 155, 207, 63 (1997). [81] W. Konigsberg and G . N. Godson, Proceedings of the National Academy of Sciences of the United States of America 5 0 , 687 (1983). [82] G. F . T . Chen and M . Inouye, Nucleic Acids 1 8 , 1465 (1990). [83] E . G . Strauss, J . H. Strauss, and A . J . Levine, in F u n d a m e n t a l Virology, edited by B. N. Fields, P. M . Knipe, and P. M . Howley (Lippincott-Raven, Philadelphia, 1996), p. 141. Bibliography 156 [84] R. Y . Tsien, Annual Review of Biochemistry 6 7 , 509 (1998). [85] A . F . Sauer-Budge, J . A . Nyamwanda, D. K . Lubensky, and D. Branton, Physical Review Letters 9 0 , 238101 (2003). [86] G . H. Vineyard, Journal of the Physical Chemistry of Solids 3, 121 (1957). [87] L. Z. Song, M . R. Hobaugh, C. Shustak, S. Cheley, H . Baley, and J . E . Gouaux, Science 2 7 4 , 1859 (1996). [88] P. G. de Gennes, Scaling Concepts Ithaca, New York, 1979). i n P o l y m e r Physics (Cornell University Press, [89] A . Hanke and R. Metzler, Journal of Physics A 3 6 , L473 (2003). [90] G . Altan-Bonnet, A . Libchaber, and O. Krichevsky, Physical Review Letters 9 0 , 138101 (2003). [91] J. Chuang, Y . Kantor, and M . Kardar, Physical Review E 6 5 , 011802 (2002). [92] O. Beckstein and S. P. Sansom, Proceedings of the National Academy of Sciences of the United States of America 12, 7063 (2003). [93] S. Iijima, Nature 354, 56 (2001). [94] Y . W. Tang, K. Y . Chan, and I. Szalai, Journal of Physical Chemistry B 108, (2004). 18204 [95] D. Nicholson and N. Quirke, Molecular Simulation 2 9 , 287 (2004). [96] A. Striolo, A . A . Chialvo, K. E . Gubbins, and P. T . Cummings, Journal of Chemical Physics 122, 234712 (2005). [97] A . Striolo, K . E . Gubbins, A . A. Chialvo, and P. T . Cummings, Molecular Physics 102, 243 (2004). [98] A . Striolo, A . A . Chialvo, P. T . Cummings, and K. E . Gubbins, Langmuir 19, 8583 (2003). [99] S. Vaitheeswaran, J . C. Rasaiah, and G . Hummer, Journal of Chemical Physics 121, 7955 (2004). [100] J. C. Liu and P. A . Monson, Langmuir 2 1 , 10219 (2005). [101] W. Y . Lo and K. Y . Chan, Molecular Physics 4 , 745 (1995). [102] M . Lee and K. Y . Chan, Chemical Physics Letters 275, 56 (1997). [103] M . Lozada-Cassou, W . Olivares, and B. Sulbaran, Physical Review E 53, 522 (1996). [104] S. H . Lee and J . C. Rasaiah, Journal of Physical Chemistry 1 0 0 , 1420 (1996). Bibliography 157 [105] M . P. Allen and D. J . Tildesley, Computer Simulation Press, Oxford, United Kingdom, 1987). of Liquids (Oxford University [106] H. J. C. Berendsen, J . R. Grigera, and T . P. Straatsma, Journal of Physical Chemistry 91, 6269 (1987). [107] J. C. Shelley, G. N. Patey, D. R. Berard, and G. M . Torrie, Journal Chemical Physics 107, 2122 (1997). [108] G. J. Tjatopoulous, D. L. Feke, and J. A. Mann, Journal of Physical Chemistry 92, 4006 (1988). [109] H. Goldstein, Classical 1980), 2nd ed. Mechanics [110] D. A . McQuarrie, Statistical 1976). (Addision-Wesley, Mechanics Reading, Massachusetts, (Harper Collins, New York, New York, [Ill] J. P. Valleau and L. K . Cohen, Journal of Chemical Physics 72, 5935 (1980). [112] I. Benjamin, Journal Chemical Physics 95, 3698 (1991). [113] K . Koga, R. D. Parra, H . Tanaka, and X . C. Zeng, Journal of Chemical Physics 113, 5037 (2000). [114] H. Tanaka and K . Koga, Journal of Chemical Physics 123, 094706 (2005). [115] G . Gao, T . Cagin, and W. A . Goddard, Nanotechnology 9, 184 (1998). [116] M . Yoon, S. Berber, and D. Tomanek, Physical Review B 71, 155406 (2005). [117] L. E . Reichl, A Modern Course in Statistical York, New York, 1998), 2nd ed. Physics (John-Wiley & Sons, New [118] G. Schonherr and G . M . Schiitz, Journal of Physics A 37, 8215 (2004). [119] G . Schonherr, Physical Review E 71, 026122 (2005). [120] P. A . Pury and M . O. Caceres, Journal of Physics A 36, 2695 (2003). 158 Appendix A A n implementation of the Bortz-Kalos-Lebowitz M o n t e Carlo Algorithm Here we present more details on the Bortz-Kalos-Lebowitz (BKL) [67] continuous time algorithm briefly described in Chapter 2. As stated in Chapter 2, the B K L algorithm for evolving the state of a stochastic system involves three steps. Given a system in a state Sk at simulation time t, the steps are 1. Compute the total transition rate (R) out of state Sk, where R — Ei=oi#fc fcJr Any state Si associated with a non-zero is accessible from state Sk in a single transition. 2. Select the transition to state Si with probability rki/R. 3. Increment the total time spent in state Sk by 1/R, record any statistical information of interest such as the current, update the system to state Si and return to step 1. One benefit of this algorithm is immediately apparent; every iteration of the algorithm is guaranteed to evolve the system to a new state. This is of critical importance for efficient simulations of the T A S E P in, for example, the high density phase where the movement of the majority of the particles is blocked, injection is unlikely, and extraction is extremely slow. Competing algorithms such as the Metropolis method [66] typically perform very poorly under these conditions as they rely on the rejection of attempted state transitions in order to determine the typical lifetime in any particular state Si. Nonetheless, while the B K L algorithm is capable of maintaining a high level of efficiency throughout the T A S E P phase diagram, there are two potentially time consuming steps in the algorithm; step one where the total rate of all possible moves is computed, and step two where a Appendix A. A n i m p l e m e n t a t i o n of t h e B o r t z - K a l o s - L e b o w i t z M o n t e C a r l o A l g o r i t h m 159 (1) s = 2 . 2 < 5 (2) 5 = 2.2 < 3 (3) 5 = 2 . 2 >0 Figure A . l : A sample T A S E P and the heap structure used in the BKL-based simulation. The value of each of the leaves equals the rate at which the associated particle move can take place. Red shaded leaves with a value of zero indicate moves that are impossible in the current state of the lattice. Listed in the column on the right are the values of s, the random selector in Alg. 6. The green hatched nodes are the nodes visited by Alg. 6 while selecting the next move, and the solid green node (Z ) is the selected move. 4 specific move is selected. Fortunately through an appropriate choice of data structure, both steps can be done efficiently. The data structure we employ in our BKL simulation is the so-called binary heap. example can be viewed in Fig. A . l . In Fig. A . l each circle represents a single node An of the heap. Each node in the heap is connected to a maximum of three other nodes. If we consider any single node in the heap n , then any nodes on a higher level in the heap x that are connected to node may are said to be children of m . have a maximum of two children which we term the left In a binary heap, a node child and the right child. Similarly, a node n may be connected to a maximum of one node that is at a lower level x in the heap, and is called the p a r e n t of m . Nodes that have no children are called of the heap, while the single node that has no parent is called the root or head leaves node. Additionally, in a binary heap every node has at least one value associated with it, and in Fig. A . l this value is represented by the number displayed inside each circle. For a Appendix A. An implementation of the Bortz-Kalos-Lebowitz Monte Carlo Algorithm 160 data structure to be classified as a heap all the values stored in each node must satisfy the heap property: the value stored in any non-leaf node rii must be greater than or equal to the values stored in the children of rij. When using a binary heap to implement the B K L algorithm for a T A S E P simulation, we specialize the heap property and require that the value stored at any non-leaf node be equal to the sum of the values of the node's children. The number of leaves in the binary heap is equal to iV + 1 where N is the number lattice sites in the TASEP. Specifically, each leaf lj stores the hopping rate of the particle located at site j. If no particle is present at site j, then lj = 0. The one extra leaf in the heap is assigned the current value of the entrance rate. Thus if a particle is at site 1 and prevents the entrance of a new particle, / try = 0, otherwise /entry = en Once the rates of all the possible particle movements have been assigned to the leaves of the heap, the heap property then guarantees that the value of the head node will equal R, the total transition rate out of the current T A S E P state. Also note that the number of levels in the binary heap is [\og (2N + 1)J + 1. 2 Having completely defined the binary heap used in the B K L method, we can now develop the algorithms for selecting the state transitions and updating the T A S E P state. First note that for a T A S E P with N sites, the associated binary heap contains 2N + 1 floating-point values. For values of N ~ 10 , it is convenient to store the entire binary 4 heap in a single array of length 2N + 1. Algorithms. 1 through 4 are then simply utility functions for locating parents and children in the heap. Algorithm 5 describes how to change a rate in an existing heap where every node already satisfies the heap property. Given a lattice site i the algorithm locates the leaf corresponding to the lattice site, and then traverses the levels of the heap from the leaf to the head node, ensuring that the heap property is satisfied. Algorithm 6 shows how, given a properly constructed heap, the next event is selected. The algorithm relies on a function R A N D O M that returns a random number s drawn with uniform probability in the range 0 to 1.0 exclusive of the end points. The uniform deviate s is then multiplied by the total rate for the T A S E P to escape its current state. Essentially the algorithm then executes a binary search for the transition with a rate r such that Efc=o fc < r { s < J2k=o k- The algorithm then returns the index r of the lattice site associated with the selected transition. Finally in Alg. 7 we outline Appendix A. An implementation of the Bortz-Kalos-Lebowitz Monte Carlo Algorithm 161 the main loop of the B K L simulation program. A n auxiliary function C O M P U T E R A T E is used that returns the movement rate of a particle located at lattice site i given the current configuration of the lattice. In the event that the lattice site is empty or the particle at the lattice site is blocked from moving, the function returns 0. As all the algorithms for manipulating the heap data structure involve traversing the heap between the leaves and the head node, these algorithms all run in 0(log (A )) time, r 2 while the storage arrays for both the heap and lattice are Q(N) in space. These highly advantageous scalings mean that this implementation of the B K L algorithm can simulate TASEPs approximately 10 sites with extremely good performance. While not absolutely 4 necessary for the standard TASEP, this heap based implementation of the B K L algorithm becomes extremely useful when more complex variants of the T A S E P are employed. Appendix A. An implementation of the Bortz-Kalos-Lebowitz Monte Carlo Algorithm Algorithm 1 LCHILD(Z) Requires: A n integer i > 1 corresponding to the array index of a node (rij) Returns: The array index of the left child of the node ri; 1: return 2i Algorithm 2 RCHILD(Z) Requires: A n integer i > 1 corresponding to the array index of a node (n*) Returns: The array index of the right child of the node rii 1: return 2i + 1 Algorithm 3 PARENT(Z) Requires: A n integer i > 1 corresponding to the array index of a node (n*) Returns: The array index of the parent of the node n; return [i/2\ 162 Appendix A. An implementation of the Bortz-Kalos-Lebowitz Monte Carlo Algorithm 163 A l g o r i t h m 4 L E A F I N D E X ( Z , N) Requires: A n integer 1 < i < N corresponding to the index of a lattice site (i) Requires: A n integer equal to the length of the T A S E P lattice (A'') Returns: The array index of the leaf corresponding to lattice site i 1: return N + i + 1 A l g o r i t h m 5 U P D A T E H E A P ( r a t e , i, Heap) Requires: A numeric value for the movement rate at lattice site, (i) Requires: The index of the lattice site, (i) Requires: A floating-point array of length 2N + 1, (Heap) Ensures: A binary heap for a T A S E P with the correct movement rate stored for lattice site i 1: 2: 3: 4: 5: 6: 7: 8: I N T E G E R j <— L E A F l N D E X ( i ) Heap[j] <— rate j <— P A R E N T ( J ' ) while j > 1 do Heap[j] <- Heap[LCmLD(j)] j <— P A R E N T ( j ) end while return + Feap(RCHlLD(j)] Algorithm 6 SELECTTRANSITION(A^, Heap) Requires: The number of sites in the T A S E P lattice, (N) Requires: A floating-point array of length 2N + 1, (Heap) Returns: The lattice site where the next event occurs l: I N T E G E R i <— 1 2: I N T E G E R k <— N + 1 {Location of the first leaf in the heap array} 3: F L O A T s <— R A N D O M ( ) x Heap{l] { R A N D O M Q generates a uniform ordinate in (0,1.0)} 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: while i < k do if s < Heap[LEFTCHILD(i)] then i <—LEFTCHILD(Z) else s <— (s - i / e a p [ L E F T C H I L D ( i ) ] ) i ^RlGHTCHILD(z) end if end while i <— (i — k) {Injection events return 0 } return i Appendix A. An implementation of the Bortz-Kalos-Lebowitz Monte Carlo Algorithm 164 A l g o r i t h m 7 S A M P L E S I M U L A T I O N ( Q ; , @,p, N, Lattice, Heap, M) Requires: A numeric value for the injection rate, (a) Requires: A numeric value for the extraction rate, (f3) Requires: A numeric value for the hopping rate, (p) Requires: Requires: Requires: Requires: An integer equal to the number of sites in the T A S E P lattice, (N) An integer array of length N with each element equal to 0, (Lattice) A floating-point array of length 2Af + 1 with each element equal to 0, (Heap) A n integer equal to the maximum number of simulation iterations, (M) INTEGER i INTEGER j FLOAT r F L O A T T {A variable to record the simulation time} T <- 0 {Initialize the B K L heap with the injection rate} U P D A T E H E A P ( a , 0, Heap) for i = 0 to M do j = S E L E C T T R A N S I T I O N ( A ' , Heap) {Select the site j where the next move occurs} {Now that we know the next move and the total escape rate, call functions to record statistical information here} T <— (T + 1/Heap[l}) {Increment the total simulation time} if 0 < j < N then Lattice[j] <— 0 Lattice[j + 1] <— 1 U P D A T E H E A P ( 0 , j, Heap) r <— C O M P U T E R A T E ( J + 1, a, (5, p, N, Lattice) U P D A T E H E A P ( 7 " , j + 1, Heap) r <— C O M P U T E R A T E ( J - 1, a, (3,p,N, Lattice) UPDATEPlEAP(r, j - I, Heap) else if j' = 0 then / U P D A T E H E A P ( 0 , 0, Heap) Lattice[l] <— 1 r <— C O M P U T E R A T E ( 1 , a, fi,p, N, Lattice) UPDATEFfEAP(r, 1, Heap) else U P D A T E H E A P ( 0 , N, Heap) Lattice[N] <- 0 r <— C O M P U T E R A T E ^ - I, a, j3,p, N, Lattice) U P D A T E H E A P ( r , N - 1, Heap) end if end for
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Transport and structure in nanoscale channels
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Transport and structure in nanoscale channels Lakatos, Gregory William 2006
pdf
Page Metadata
Item Metadata
Title | Transport and structure in nanoscale channels |
Creator |
Lakatos, Gregory William |
Date Issued | 2006 |
Description | Driven by the rapidly advancing fields of nano- and biotechnology, there has been an explosion of interest in molecular transport and structure formation on small length scales. A canonical model for the transport of particles along one dimensional pathways in nanoscale channels is the Totally Asymmetric Simple Exclusion Process (TASEP). After introducing the standard TASEP, modifications of the TASEP designed to increase its utility in modeling biological transport processes are described. One variant of the TASEP is particularly suitable for modeling protein translation, and the results of using this variant to investigate the effects of slow-codons on the translation process are discussed. A related topic is the voltage-driven translocation of DNA hairpins through membrane-embedded nanopores. Motivated by recent experiments, a stochastic model is developed that couples the translocation and dehybridization of the DNA hairpin. This model is used to explore the behaviour of the mean translocation time of hairpins as a function of driving voltage, and two translocation mechanisms are identified and discussed. Finally, the adsorption and equilibrium structures of water in the interior of ion-bearing nanoscale pores are considered. The behaviour of water and ions under confinement is critical to the functioning of biological ion channels and nanoporous filters. Here, the adsorption isotherms of water are examined, and the layered structures formed by the confined water are described. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-01-17 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0085456 |
URI | http://hdl.handle.net/2429/18543 |
Degree |
Doctor of Philosophy - PhD |
Program |
Physics |
Affiliation |
Science, Faculty of Physics and Astronomy, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2006-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_2006-200070.pdf [ 20.2MB ]
- Metadata
- JSON: 831-1.0085456.json
- JSON-LD: 831-1.0085456-ld.json
- RDF/XML (Pretty): 831-1.0085456-rdf.xml
- RDF/JSON: 831-1.0085456-rdf.json
- Turtle: 831-1.0085456-turtle.txt
- N-Triples: 831-1.0085456-rdf-ntriples.txt
- Original Record: 831-1.0085456-source.json
- Full Text
- 831-1.0085456-fulltext.txt
- Citation
- 831-1.0085456.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
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"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0085456/manifest