UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Transport and structure in nanoscale channels Lakatos, Gregory William 2006

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

Item Metadata

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

Transport and Structure in Nanoscale Channels by Gregory W i l l i a m Lakatos A T H E S I S S U B M I T T E D I N P A R T I A L F U L F I L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F D O C T O R O F P H I L O S O P H Y in T H E F A C U L T Y O F G R A D U A T E S T U D I E S (Physics) T H E U N I V E R S I T Y O F B R I T I S H 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 explo-sion 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 Total ly Asymmetric Simple Exclusion Process ( T A S E P ) . After introduc-ing 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 par-ticularly 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 membrane-embedded nanopores. Motivated by recent experiments, a stochastic model is developed that couples the translocation and dehybridization of the D N A 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. i i i Contents Abstract i i Contents i i i List of Tables v List of Figures v i Preface x iv Acknowledgements xvi i 1 Introduction 1 2 The Totally Asymmetric Simple Exclusion Process 4 3 Totally Asymmetric Exclusion Processes with Finite Size Particles . . 12 3.1 Background 12 3.2 Mean Fie ld Theories 13 3.2.1 Simple Mean Fie ld Methods 13 3.2.2 Refined Mean Fie ld Theory for the M a x i m a l Current Phase . . . . 17 3.2.3 Refined Mean Fie ld Theory for Boundary-l imited Currents 19 3.2.4 Matching Mean Fie ld Phase Boundaries 22 3.3 Monte Carlo Simulations 23 3.3.1 Currents 24 3.3.2 Particle Densities . 27 3.3.3 Phase Boundaries 29 3.4 Summary 33 4 A Totally Asymmetric Exclusion Process with Periodic Structure . . . 34 4.1 Background 34 4.2 Mode l and Methods 35 4.3 Mean F ie ld Theories 36 4.3.1 Simple Mean Fie ld Methods 36 4.3.2 Refined Mean Fie ld Methods 38 4.4 Monte Carlo Simulations 48 4.4.1 Currents 48 4.4.2 Densities 52 4.4.3 Phase Diagram . 56 Contents iv 4.5 Other Periodic Arrangements 56 4.6 Summary 61 5 Clustered Bottlenecks in m R N A Translation and Protein Synthesis . . 63 5.1 Background 63 5.2 Mode l 65 5.3 Results and Discussion 69 5.4 Summary 74 6 First Passage Times of Driven D N A Hairpin Unzipping 76 6.1 Background 76 6.2 Model 78 6.3 Results and Discussion 86 6.4 Summary • • 95 7 Water Adsorption in Ion bearing Nanopores 97 7.1 Background 97 7.2 Mode l and Simulation Methodology 98 7.2.1 Molecular Models 98 7.2.2 Simulation Methodology 104 7.3 Results and Discussion 110 7.3.1 Adsorpt ion Isotherms 110 7.3.2 F i l led State Structures 114 7.3.3 Adsorption In Ion-Bearing Pores 116 7.3.4 Pressures 132 7.3.5 Grand Potentials : 134 7.3.6 Boundary Effects 139 7.4 Summary 143 8 Conclusions and Suggestions for Future W o r k 146 Bibliography 152 A A n implementation of the Bortz-Kalos-Lebowitz Monte Carlo Algo-rithm 158 List of Tables 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 79 6.2 List of fitting parameters used in the two dimensional stochastic model. . 85 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. . 87 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 F ig . 6.3(a) 88 7.1 The 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]. The carbon param-eters were used for the wall potential given in E q . 7.2 99 vi List of Figures 2.1 (a) The 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 wi th an attempt rate (3. No more than one particle can occupy a lattice site at a time, (b) The 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 : The density profile along the T A S E P lattice. In the entry and exit l imited phases two density profiles are possible depending on the values of (a,/3) relative to the a = 1 — P line (the dashed line in the plot) 6 3.1 Entry and exit mechanisms, (a) Par t ia l entry (b) complete entry into the first d sites (c) partial exit (d) complete exit as the particle reaches the last site 15 3.2 (a) A typical configuration of particles associated wi th 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 with each particle is shown with a blue rectangle (b) A typical configuration for the calculation of Z(n, L — 1). 17 3.3 State enumerations for particles of size d = 4 near the lattice boundaries, (a) The d+l distinct states at the entrance region, (b) The 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 in the states displayed in (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 wi thin d moves of state Si are included 19 3.4 The currents wi thin 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) The entry limited current JL as a function of d for (3 = 10 and a = 0.1. (b) The maximal current J m a x for a = ft = 10. (c). The exit l imited current JR for a = 10 and (3 = 0.1 25 3.5 Currents as a function of a or j3 for d = 2 ,4, and 8. (a) The entry limited to maximal current profile with (3 = 10. (b) The exit l imited to maximal current profile wi th a = 10. The 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) 26 3.6 Boundary densities in each of the three current regimes. The solid lines show the predicted densities from the refined mean field approximations (Eq. 3.24) 28 Lis t of Figures v i i 3.7 Representative scaled density (aid) profiles in the three current regimes, (a) The entry rate l imited regime (a = 0.1, (3 = 1). (b) The maximal current regime (a = (3 = 1). (c) The exit rate l imited regime (a = l , /3 = 0.1). The boundary regions are shown in the insets. The strong splitting of the densities into high density and low density branches, is a product of the close-packing of the particles in the exit l imited regime, and the fact that particle positions are determined by the location of the particles' left edges. 30 3.8 The phase diagram for the T A S E P with d = 3 particles derived from the refined mean field results (Eq. 3.24) and the simple mean field results (Eq. 3.8). The points are the transitions estimated from Monte Carlo simulations. 31 3.9 (a) The discontinuity in the second derivative of the steady-state current wi th respect to a or B at the transition between the maximal current and the entry limited regimes. The simple mean field approach (Eq. 3.9) pre-dicts an increasing discontinuity with d, while the refined mean field ap-proach (Eq. 3.25) predicts a jump that decreases wi th d. (b) The critical values a* and B* at which the boundary limited to maximal current tran-sitions occur (Eqs. 3.8 and 3.24) 32 4.1 The periodic dual-rate totally asymmetric exclusion model. A t every T lattice sites, the particle movement rate is p2. 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 35 4.2 The 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 . 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 « e f f = P2&4 on the left, and an effective extraction rate of Be$ = ^2(1 — en) on the right 43 4.3 The steps in 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 wi th the corresponding decimal value (i.e. lattice occupancy O i l —» state 3). The 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 in 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 in the figure). Cal l ing 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 ina l ly the transitions between each 0-state and each 1-state generated by injection at the left edge of the lattice (dotted arrows) are enumerated. . . 44 List of Figures vi i i 4.4 (a) The 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 maximal current phase (a = 10, /? = 10). (b) Max imal currents for dual-rate T A S E P s of various periods. The F S M F estimates shown by solid lines were produced using a single period of the lattice as the finite segment. Whi le 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 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 with increasing T can be partially reduced by increasing the segment length 49 4.5 (a) Current profiles along the a direction in the T — 2 dual-rate T A S E P phase plane (B — 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 T A S E P s of various periods wi th p2 = 0.25. The F S M F results (lines) were produced using an L = T F S M F approach (/? = 10) 51 4.6 Monte Carlo density profiles for a T = 2 T A S E P w i t h p 2 = 0.1. Comparison wi th F ig . 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 l imit ing boundary, consistent wi th the assumptions leading to Eqs. 4.18 and 4.19 52 4.7 Max ima l 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 assump-tion produces the largest errors at small values of p 2 where the density correlation between adjacent sites is expected to be large 53 4.8 The bulk density profiles for the T = 5 and T = 9 dual-rate T A S E P s wi th p2 = 0.2, in the maximal 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 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 54 List of Figures ix 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 E q . 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 E q . 4.20. Al though the solution for cr2 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 o2 profiles predicted by the two mean field methods at small values of p2 55 4.10 (a) The phase diagram for the T = 2 two-rate T A S E P . 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 two-rate T A S E P at various periodicities, wi th p2 = 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 57 4.11 Variations on the periodic dual-rate T A S E P . The 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 (p2 = 0.25). The lines are the best-fit phase boundaries for T A S E P s with values of N that satisfy Eq . 4.28. The points are the Monte Carlo predictions for the phase boundaries of T A S E P s wi th values of N that violate Eq . 4.28. Note that in both 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 /?*. Addit ionally, 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, and share a common value of (3* 60 5.1 The m R N A translation/protein synthesis process. Ribosomes move unidi-rectionally along the m R N A as t R N A s (not shown) deliver the appropri-ate amino acid to the growing protein. Codons wi th 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 m R N A translation. Slow-codons are represented by lattice sites with a reduced particle movement rate p2 < p\. 64 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. . 66 L i s t of F i g u r e s x 5 . 3 (a) Comparison of steady-state currents J\{p2) for a single defect derived from M C simulations wi th those obtained from a L-segment F S M F method. For L > 5 , the F S M F results are within 2 % of those from M C simula-tions. The boundary injection and extraction rates were not rate l imiting (a = 8 = 1 0 ) . (b) Further reduction of the steady-state current as succes-sive, identical defects are added. The first few defects cause most of the reduction in current. Note that the ribosome current varies significantly as P2 is changed. Thus for clarity, the currents in the figure have been normalized by the slow-codon hopping rate ( j ^ ) - Also note that for small P2, J2/P2 ~ 1 / 2 consistent with Eq . 5 . 3 6 8 5 . 4 (a) Steady-state currents across a lattice wi th 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 2 = 0 . 1 5 ) of various spacings k from F S M F computations. The circled points corre-spond to the defect positions for k = 6 . For larger k, density boundary layers decay, allowing particle build-up behind the second barrier, enhanc-ing the current 7 2 5 . 5 (a) The dependence of the normalized steady-state current J2(P2, k)/Ji{.P2) as a function of the defect separation, (b). The upper bound for the mean current of a lattice wi th m defects randomly distributed within the central N = 3 0 0 sites. A lower bound as m — > N is J ( m — > N)RAN —* P2 /4 7 3 6 . 1 (a) Divis ion of the a-hemolysin channel into the cis, vestibule, transmem-brane, 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 NV = 1 6 , while the number of D N A bases than can fit into the transmembrane region is NL = 1 2 . (b) A n example of the D N A hairpin translocation experiment. The process of D N A translocation is described by two coordinates n 0 and ri\. n0 gives the position of the first hairpin base relative to the mouth of the transmembrane part of the channel, while index rii gives the number of hairpin base pairs that have been dehybridized. . . 7 9 6 . 2 The enumeration scheme for D N A base pairs in the hairpin and the asso-ciated state space for the 2 D stochastic model. In this example N = 7 , NA — 4 , and there is a defect at base 4 . States wi th n\ = 3 are inac-cessible 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 nx have attempt frequency 00. In all cases, transitions between states wi th different values of n 0 have attempt frequency \x. States with n 0 = 2A f + A / loop + NL (a D N A strand pulled completely into the trans chamber), and n 0 = — NA (a D N A strand pulled completely into the cis chamber) are absorbing 8 0 List of Figures XI 6.3 (a) The best fit predictions of the stochastic model when both 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 primari ly by sliding back into the cis chamber. A t 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 po ly-A tail , (b) The probability of the D N A strands escaping into the cis chamber, (c) The 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 defect-bearing 2-3 free energy surface 90 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 fc^T is taken to equal the physical radius of the pore 102 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. A t the density of l iquid 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 I l l 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 (m 1.0g/cm 3 at 298K), while the vertical line gives the chemical potential of saturated water vapour at 298K ( « - 5 8 . 5 k J / m o l ) 113 7.4 Characteristic water structures in pores wi th 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 wi th an average of 1.8 hydrogen bonds per water molecule. \x\v = —56.0kJ/mol 115 7.5 Water Oxygen - Water Oxygen radial distributions functions for pores con-taining pure water, and pores containing single ions (fiw = —56.0kJ/mol). The ions produce deviations from the pure water structure that are most pronounced in pores wi th smaller radii, and increase wi th the size of the ion's charge. In the 0.74nm pore, the perturbation of the water structure resulting from the ion is small 117 7.6 F i l l i ng in the 0.44nm pore with a C l ~ ion. (a) fiw = —66kJ/mol (nw = 4) (b) fj,w = — 64kJ /mol (nw = 6) (c) the low density state at nw = —62kJ/mol (nw = 7) (d) the high density state at jiw = —62kJ/mol (nw = 11) (e) fiw = — 59kJ /mol (nw = 12) (f) fiw = —56kJ/mol (nw = 13). A t fiw = —56kJ/mol states with 12 water molecules similar to the state shown in (e) and states with 13 water molecules are equally probable. 118 List of Figures x i i 7.7 F i l l ing in the 0.44nm pore wi th a K + ion. (a) fiw = —66kJ/mol (nw = 5) (b) fjLW = - 6 4 . 0 k J / m o l (nw = 7)(c) fiw = - 6 2 . 0 k J / m o l (nw = 9) (d) fiW = - 6 1 . 0 k J / m o l (nw = 10) (e) fiw = - 5 9 . 0 k J / m o l (nw = 13) (f) HW = - 5 6 k J / m o l (nw = 14) 119 7.8 F i l l ing in the 0.54nm pore with 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 (nw = 31) (e) Ltw = - 5 9 k J / m o l (nw = 34) (f) fiw = - 5 6 k J / m o l (nw = 40). 120 7.9 F i l l ing in the 0.54nm pore with 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) LLW = —58.9kJ/mol (nw = 31) (f) JJLW = —56kJ/mol (nw = 37) 121 7.10 F i l l ing in the 0.64nm pore with a C l ~ ion. (a) Liw = —66kJ/mol (nw = 6) (b) jjw = —64kJ/mol (nw = 7) (c) the low density state at /J,W — —60.6kJ/mol (nw = 19) (d) the high density state at fj,w = —60.6kJ/mol (nw = 50) (e) fiw = —59.8kJ/mol (nw = 53) (f) fiw = —56kJ/mol (nw = 63) 122 7.11 F i l l ing in the 0.64nm pore with 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) fiw = —58.9kJ/mol (nw = 52) (f) fiw = —56kJ/mol lnw = 60) 123 7.12 F i l l ing in the 0.74nm pore with a C l ~ ion. (a) fiw = —66kJ/mol (nw — 5) (b) fiw = - 6 4 k J / m o l (nw = 6) (c) fiw = - 6 0 . 3 k J / m o l (nw = 20) (d) immediately prior to filling fiw = —59.9kJ/mol (nw — 25) (e) immediately after filling fiw — —59.8kJ/mol (nw = 83) (f) fiw = —56kJ/mol (nw = 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. 124 7.13 F i l l ing in the 0.74nm tube with a K + ion. (a) fiw — —66kJ/mol (nw = 6) (b) fiw = —64kJ/mol (nw = 7)(c) immediately prior to filling fiw — —59.7kJ/mol (nw = 24) (d) immediately after filling fiw = —59.4kJ/mol (nw = 82) (e) fiw = —57.5kJ/mol (nw = 88) (f) fiw = —56kJ/mol (nw = 92). In (f) a second inner layer of water is beginning to form 125 7.14 The total water-ion interaction energy (closed points), and the total particle-wall interaction energy (open points) for the monovalent ions in (a) 0.44nm, (b) 0.54nm, (c) 0.64nm, and (d) 0.74nm pores 126 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 126 7.16 Radia l probability density profiles for C l ~ pores wi th 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 127 List of Figures xi i i 7.17 Radia l probability density profiles for K + pores wi th radii (a) 0.44nm (b) 0.54nm (c) 0.64nm (d) 0.74nm. Note that the filled state hydrogen distri-butions exhibit increased density at large radii compared to the profiles of F ig . 7.16 128 7.18 Radia l distribution functions for K + and C P in pores wi th radii equal to 0.44nm and 0.54nm. In both cases \i\v — —59.0kJ/mol 129 7.19 (a) The radial pressures in K + (closed symbols) and C P (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 . . 133 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 PL over the range of displayed values of L is much smaller than the variation in the K + bearing pore ([xw = —57.0kJ/mol) 136 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 138 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 ( k J / m o l • K ) of particles 139 7.23 (a) The adsorption isotherms for K + bearing pores wi th radius R = 0.54nm with different lengths. The (b) ion, (c) water oxygen, and (d) water hy-drogen radial probability profiles. The points are for the L — 1.5nm pore while the lines are for the L = 3.1nm pore (LIW — —58.2 kJ /mol ) 140 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 wi th the same Lennard-Jones parameters as K + . In all cases the pore radius was R = 0.28nm and the pore was filled wi th n w = 34 water molecules 141 A . l A sample T A S E P and the heap structure used in the B K L - b a s e d simula-tion. The value of each of the leaves equals the rate at which the associated particle move can take place. Red shaded leaves wi th a value of zero indi-cate 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 A l g . 6. The green hatched nodes are the nodes visited by A l g . 6 while selecting the next move, and the solid green node (Z4) is the selected move 159 xiv Preface The content of this thesis has been based, in part, on the following articles: 1. "Totally asymmetric exclusion processes with particles of arbitrary size", G . Lakatos, T . Chou, J. Phys. A 36 (2003) 2027 (Chapter 3) The original idea for the research presented in this paper was due to Professor Tom Chou ( U C L A ) . The design of the research program was suggested by Professor Chou, while the methodologies and results were developed by G . Lakatos wi th guidance from Professor Chou. The final manuscript was written by G . Lakatos with revisions and additions by Professor Chou. 2. "Steady-state properties of a totally asymmetric exclusion process with periodic structure", G . Lakatos, T . Chou, and A . Kolomeisky Phys. Rev. E 71 (2005) 011103 (Chapter 4) The idea for the research presented in this paper was originated by Professor Anatoly Kolomeisky (Rice) and Professor Tom Chou ( U C L A ) . The design of the research program, the methodologies, and the presented results were developed by G . Lakatos with guidance and suggestions from Professors Chou and Kolomeisky. The final manuscript was written by G . Lakatos with revisions and additions from Professors Chou and Kolomeisky. 3. "Clustered Bottlenecks in m R N A Translation and Protein Synthesis" ,T. Chou, G . Lakatos Phys. Rev. L e t t . 93 (2004) 198101 (Chapter 5) The idea for the research presented in this paper is due to Professor Tom Chou. Professor Chou 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 wi th guidance and suggestions from Professor Chou. Professor Chou co-wrote the text wi th 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) The research was proposed by G . Lakatos and Professors Bergersen and Patey. The design of the research program was developed by G . Lakatos and Professor Patey. The research and results are due to G . Lakatos. The manuscript was written by G . Lakatos wi th suggestions from Professors Patey, Chou, and Bergersen. 5. "Water Adsorption in Ion Bearing Nanopores", G . Lakatos, G . N . Patey, To be submitted to the Journal of Chemical Physics (Chapter 7) The idea for the research presented in this paper was conceived by G . Lakatos and Professor Patey. The design of the research program was developed by G . Lakatos and Professor Patey. The research and results are due to G . Lakatos with guidance and suggestions from Professor Patey. The manuscript was written by G . Lakatos with suggestions from Professor Patey. As part of the body of software used in the preparation of this thesis, the following third-party numerical analysis codes were employed 1. The Arno ld i matrix diagonalization routines [1] employed in Chapters 4, 5, and 6. The implementation used is included wi th the commercial software package M A T L A B [2]. 2. The simplex multi-dimensional minimization scheme [3] employed in Chapter 6. The implementation used is included with the commercial software package M A T L A B . 3. The Mul t ip le-Fold ( M F O L D ) 2-state hybridization server used in Chapter 6 [4]. Preface xvi A l l other simulation and analysis codes, including the Monte Carlo simulation programs used in Chapters 3, 4, 5 and 7, were developed by G . Lakatos. This includes extensive code and scripting developed to drive the third-party numerical analysis codes enumerated above. XVII Acknowledgements Writ ing 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. They 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 in particular Professor Tom Chou, who has provided helpful suggestions and good advice. The University itself has provided a wonderful environment in 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 in molecular transport pro-cesses 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, sev-eral 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 al l l iving organisms. Examples include the movement of so-called motor proteins such as myosin [5] that is ultimately responsible for all muscular movement in multi-celled life. Myosin molecules, along with many other motor proteins, have been found to move in one direction along long polymer-like 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 expres-sion, and protein synthesis [7]. Here we find multiple distinct transport processes that are critical elements in the mechanism that converts the information stored in an organism's D N A into functional proteins. As 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]. This 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. The movements of motor proteins and ribosomes have been investigated using a va-riety 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 adenosine-triphosphate molecule [9], as well as the binding of amino-acid carrier molecules ( tRNAs) to ribosomes. Similarly, molecular dynamics simulations [10] have been used to study the conformation changes involved in 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 trans-port systems, and are difficult to generalize to other biological transport processes. This latter point is particularly acute given that many biological transport processes share a number of pertinent features. For example, the movement of both ribosomes and myosin molecules occurs along stand-like filaments and so can be treated as one-dimensional. Additionally, in both cases the molecules only move in one direction along the filaments. Finally, myosin molecules and ribosomes both move in 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 in a variety of contexts ranging from the study of highway traffic [19, 20, 29], to the foraging behaviour of ants [35]. As the T A S E P model forms the basis of the first part of this thesis, Chap-ter 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 util i ty in modeling biological transport processes. The extensions of the T A S E P and the solution techniques developed in Chapters 3 and 4, are then applied in 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 with its position along the length of the m R N A . Another example of transport in a biologically derived system is the translocation of D N A and R N A through membrane embedded nanopores. The translocation of D N A and R N A through membrane pores occurs naturally in all organisms containing a cell nucleus, and is a crit ical 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 in the development of several models of the thermodynamics of the translocation process. This 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 . Combining 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. The model is used to investigate the variations in 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 in Chapter 7, where a study of water adsorption in narrow cylindrical pores is presented. The 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 in 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 in the ion mobility by a factor of two to four compared to the bulk mobilities [60]. Addit ionally, there has been significant recent interest in util izing 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 Monte 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, in 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 wi th open boundaries at either end (see F ig . 2.1(a)). Particles exactly one lattice spacing in size can enter the lattice at the left end with 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 wi th 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 with an attempt rate B. In the variant of the T A S E P discussed in 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. This 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]. As 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 both the transported objects (the particles) and the steps they take, and it's unidirectional as all the particles move left to right. As wil l 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 in the lattice per unit time, and ct j , the particle density at any particular lattice site i. The density at a site i was defined to be the probability of finding the site occupied C h a p t e r 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 by a particle. 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. Whi le the early investigations of the model did not produce an exact solution, the T A S E P was shown to exist in one of three solution regimes, or phases, depending on the relative values of the entry and exit rates (see F ig . 2.1(b)). In the case where a < ft, and a < p/2 the T A S E P was in a low density entry-limited regime. In this case, getting particles into the lattice was the step that l imited the overall rate of particle transport. Similarly, when P < a and P < p/2, the T A S E P was in a high density exit-limited regime where removing particles at the right end of the lattice limited the current. F ina l ly when a and P were both greater than p/2 the T A S E P was in 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. The boundaries between the phases shown in Fig . 2.1(b) are sharply marked by discontinuities in measurable properties of the T A S E P . Specifically, at the transition between the entry and exit l imited 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 all show a jump discontinuity. Similarly, at the boundary limited to maximal 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 in the current and thus the entry to exit l imited 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 in the study of the T A S E P was the development of exact solutions by Derr ida et al. [15] and Schiitz and Domany [16]. The matrix product solution of Derrida in particular provides concise expressions for the mean densities, density correlations, and currents in the steady-state of the T A S E P . For example, the current in a Af-site T A S E P is given by the expression j = RN-I(p/P) - RN-iip/a) ~ P RN(p/P)-RN(p/a) (2-1) C h a p t e r 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 0 T T T T T i (b) 1 J ! V 1 1 Exit Limited Phase 0.6 Figure 2.1: (a) The standard totally asymmetric simple exclusion process ( T A S E P ) . Par-ticles 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 with an attempt rate B. No more than one particle can occupy a lattice site at a time, (b) The 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 : The 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 in the plot). Chapter 2. The Totally Asymmetric Simple Exclusion Process 7 where, R M ) - " T U - 1 ) { 2 N - i ) l x < (2 2) HN(x)-2^ m , N + 1 _ j ) l * - (2-2) In the l imit of large JV, the exact expression for the current takes on one of the three phases shown in F ig . 2.1(b) and shows the same discontinuities at the phase boundaries as were known to MacDonald et al. [11, 12]. Addit ionally, 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 in the three solution regimes (Entry) a < p/2 JL = a(l - a/p), aL = a/p a < P (2-3) (Exit) 0<p/2 J R = p(l-p/p),oR = l-p/p 8 < a (Maximal) a > a* = p/2 J m a x = p /4 , a m a x = 1/2 The expressions for the current in Eq . 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- We then set OL equal to the ratio of the injection rate to the internal stepping rate, ai — ct/p. The particle current into the first site in the lattice is then JL = a(l — a/p). A similar argument gives JR. = P(l — P/p) in the exit-limited phase. In the maximal current phase we assume that the density is constant and the probability of any particular site in the lattice being occupied is uncorrelated with 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). Maximiz ing this expression gives a = 1/2 and J = p/A. As we ignore correlations between lattice sites when producing these results, these simple calculations are so-called mean field calculations. C h a p t e r 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 8 Another important feature of the T A S E P that manifests itself in E q . 2.3 is the sym-metry between the entry- and exit-limited phases. Consider the entry-limited solution in Eq . 2.3 and replace a wi th (3. This immediately returns the correct current for the exit limited phase. Similarly, the exit-limited density equals 1 — o_ , under the substitution a —> (3. This 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 in 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. The mean occupation of an interior lattice site and the first lattice site are then given by d ( x i ) = p ( ( X j - l ) - (Xi) + {XiXi+l) - ( X i - i X i ) ) (2.4) = a ( l - E i ) - p ( x i ( l - x 2 ) > . Here, we can note that the symmetry-exchange operation interchanges a and f3 and re-places Xi wi th (1 — xN-.i+i) in E q . 2.4, yielding = P ( ( x N - i ) ~ ( Z j v - i + l ) - {xN-iXN-i+i) + (xN_i+iXN_i+2)) dt , (2.5) d t = -p(xN)+p{xN-i(l-xN)). These are the correct equations for an interior lattice site x j v - 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. This 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 C h a p t e r 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 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 in 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 in the T A S E P lattice decay to zero far from the lattice boundaries. If this were true, then in the interior of the lattice mean field theories would be exact. As 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 in the region where mean field theory is exact wi l l 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 in the T A S E P are found to decay to zero as the inverse of the distance from the boundaries in the maximal current phase, and exponentially in the boundary limited phases [15]. The 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 maximal current phase, and exponentially in the boundary limited phases [15]. As mentioned above, while mean field approaches have been useful in studying the standard T A S E P , they do not successfully capture all the interesting aspects of the model. This problem is particularly acute when dealing with 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 Monte Carlo method is a means of accurately calculating the probability of finding the system in any particular state St. From a knowledge of the state probabilities, any property of the system can in 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^. Each 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 with a single index k, and then we have the following equation for the time evolution of the probability distribution for the lattice state In Eq . 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 in the interior of the lattice, then Tik = p . Equation 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. Many different Monte Carlo algorithms for determining the steady-state probabilities of the lattice states Si are available, and in 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 in 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^krki- A n y state Si associated wi th a non-zero is accessible from state Sk in a single particle movement. 2. Select and perform the transition to state Si wi th 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 in Appendix 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 wi l l be well approximated by the ratio of the time spent in state Sk to the total simulation time. Note that in 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 wi l l evolve the lattice from one state to the ,N (2.6) Chapter 2. The Totally Asymmetric Simple Exclusion Process 11 other. 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 Eq . 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. Divid ing the total number of particles that exit the lattice by the total simulation time yields the Monte Carlo 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. The 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 Monte Carlo simulations are used to investigate the steady state currents and densities in two modified T A S E P models. Chapter 3 removes the size l imitat ion 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 MacDonald '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 in 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 Particles1 3.1 Background As described in Chapter 2, the standard T A S E P consists of particles that make steps with a length equal to the particle diameter. However, in 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 in 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 with 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 with d > 1. Our approach is guided by the results of Monte 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 in 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 maximal 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 Mean 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 with N ^> d sites. The particles interact through hard core repulsion. The entry rate a, and exit rate (3 are all normalized by the internal hopping rate p. As a result, all current results wi l l 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. The 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 - O i + 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 C h a p t e r 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 P a r t i c l e s 14 with sizes greater than a single lattice site; namely that when a particle is at site i, the next (d — 1) sites are empty with 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)d. (3.2) The maximum current, and its associated density is then found by setting d J / d a = 0. This yields dd 1 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 end. Figure 3.1 shows two possibilities for the rules of entry and exit when d > 1. For example, a particle may only enter incrementally, one lattice site at a time (Fig. 3.1(a)), or it may enter completely in one step when the first d sites of the lattice are empty (Fig. 3.1(b)). However, in the large l imit , 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 in 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 in the l imit of large N. For the sake of simplicity, we wi l l assume complete entry (Fig. 3.1(b)) in all of our subsequent work. In the simple mean field method we wil l consider complete exit (Fig. 3.1(d)), while in the refined mean field method described later partial exit (Fig. 3.1(c)) wi l l be used. To find the entry-limited current in the simple mean field approximation, we enforce current continuity at the left edge of the lattice giving d d aY[(l-aj) = a1]l(l-a1+j). (3.4) 3 = 1 j = l Setting Oi = o~\ = a for all i, the occupation probability near the left entry site is then C h a p t e r 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 P a r t i c l e s 15 ( a ) ( c ) ! .... .• I I ! i 1 i 1 ! ! I ! ; 1 '• i i i -Y i i ( b ) ( d ) i r r - l r—r—r -1 i - r n r i r - m m n m ~ a P Figure 3.1: Entry and exit mechanisms, (a) Part ial entry (b) complete entry into the first d sites (c) partial exit (d) complete exit as the particle reaches the last site. o = oL = a and the entry-limited steady-state current becomes JL = a(l-a)d. (3.5) The T A S E P wi l l be in 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 JL < 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. When the current is exit rate limited we again enforce current continuity yielding: d - l d P&N-d+l J _ ( l - &N-d+l+j) = ON-d Jj£(l - CTN-d+j)- (3.6) i=l j=i As before we set ctj = a for all i. This gives JR = Pd(l-(3). (3.7) Setting Eq . 3.7 equal to Eq . 3.3 gives 0* = d/(d+ 1). The steady-state currents and densities produced by the simple mean field approach can be summarized as C h a p t e r 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 P a r t i c l e s 16 (Entry) a < - ^ - j - JL = a(l - a)d oL = a a(l-a)d <Bd(l- B) (Exit) ( 5 < ^ l JR = Pd(l-P) aR = l-P Pd{l - P) < a(l - a)d (3-8) 1 dd 1 (Maximal) a > a* = — j _ = ^ = _ _ d + 1" The currents, boundary densities, and interior densities in 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 da2, d2r da2 = A y=a(*)- da2 v=rv(*) d2r dP\ dd~ /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 wil l be compared with more accurate mean field theories and numerical simulations. C h a p t e r 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 P a r t i c l e s 17 Figure 3.2: (a) A typical configuration of particles associated with the calculation of Z(n, L) given by Eq . 3.11. The left edge of the particles is noted by a black dot, while the exclusion zone associated wi th each particle is shown with 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 wi l l 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 with locally uniform particle distributions is to consider a means of explicitly including the extended exclusion due to the length of the particles. Again 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 ig . 3.2). Here we take L » d and L <C N. The 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{Tl+l > J + d\Ti = j) = Z{n,L- 1) (3.10) Z(n, L) ' In Eq . 3.10, n is the position of the left edge of particle i and Z(n,L) is the number of ways of arranging n particles of size d in a lattice of L sites. Note that Z(n, L) is C h a p t e r 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 P a r t i c l e s 18 exactly the canonical parti t ion 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 Z(n,L- 1) P ( r l + 1 > j + d\n = j) = Z(n,L) L — (d — l ) n — n 1 — ad (3.12) L-(d-l)n l-a(d-l)' In Eq . 3.12 we have assumed the number density of particles in the lattice interior is approximately constant and set a = n/L. The steady-state current is then J = P(n = j)P(n+1 > 3 + d\n = j) = <^1^~^d_iy (3-13) The maximal current J m a x and its associated bulk density a m a x is found by maximizing J with respect to a yielding 1 an Jry Vd(Vd+l): x (3-14) (Vd+1)2' As expected, E q . 3.14 reduces to the standard T A S E P result [13, 15] when d — I. Also in the l imit of large d, J m a x and a m a x ~ 1/d. Note that a m a x is the probability of finding the left edge of a particle at any particular lattice site in the interior of the lattice. However, any particular lattice site i wi l l 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 damax. The steady-state current wi l l be the maximal 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 l imit ing steps. A s wi l l be demonstrated later in the chapter, the current and density derived from this refined mean field approach C h a p t e r 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 P a r t i c l e s 19 (a) a J i = i sl 0 0 0 0 h 1 0 0 0 S3 0 1 0 0 0 S4 0 0 1 0 0 0 S5 0 0 0 1 0 0 0 (b) i = N - 4 N - 3 N - 2 N - l N J _ i 0 0 0 1 sl I P 0 0 0 S2 ! 1 0 0 0 S3 i 0 1 0 0 S4 1 0 0 1 0 S5 p 0 Figure 3.3: State enumerations for particles of size d = 4 near the lattice boundaries, (a) The d+1 distinct states at the entrance region, (b) The 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 in (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 within d moves of state s\ are included. match extremely well wi th 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 ig . 3.3(a). As 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 C h a p t e r 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 P a r t i c l e s 20 JL = aP{n >d) = Pin = d,r2> 2d). (3.15) 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 within 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 s5 for d — 4 particles shown in F ig . 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~x + (1 — 9)d is the sum of the probabilities of all possible states of the first d lattice sites (Fig. 3.3(a)). The 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 ( r i = d , T 2 > 2 d ) = ^ ( i _ ^ - ^ i _ f l ) j . (3.17) We have assumed that 9 is approximately constant over the first d sites. The mean probability of occupation of the first d sites is then 9(l-9)d~l 9 ° L ~~ de(i - ey-1 + (i - ey ~ i + (d-i)e' ( 3 ' 1 8 ) Using E q . 3.15 we find that 9 = a, thus in the entry l imited phase a Jr, = 1 + (d- l)a a(l — a) (3.19) 1 + (d-l)a In Eq . 3.19, oi is the mean probability of finding site i occupied by the left edge of a C h a p t e r 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 P a r t i c l e s 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 ioL. For 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 wi l l split into two classes. Lattice sites with index i = N — kd for integer k > 0, are the sites where the left edge of the particles spend the majority of the time in the close packed state. As 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 ig . 3.3(b), we then apply current continuity to transitions from state s 5 to state S\. This gives the following relation for the steady-state current JR = (3P(rn = N) = P{rn = 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 — 1 is l " ' (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. The quantity ( l _ 0 ) d - i + (__ i ) 0 ( i -e)d~2, (3.22) is then the parti t ion function for the last set of low-occupancy sites N — d+l,N — d + 2,...,N-1. Solving E q . (3.20), we find (36N/{1 -6N) = 9/(1 + (d-2)9). To find another relationship between 9 and 9N, we can equate the current into and out of the state s 3 which gives 6^(1 - 9N)(1 - O ) ^ 1 = 9(1 - 9)d~2(l - 9N). Solving these two equations we find C h a p t e r 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 P a r t i c l e s 22 OR JR = 1-/3 l + {d 0 ( 1 --l)p' P) l + {d-l)P' (3.23) and mean occupations at the low density sites of r rw-i = P&R- Again the values of o here are probabilities of finding the left edge of a particle in 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 Matching Mean Field Phase Boundaries The different approximations for the steady-state currents and densities cover all allowable values of the parameters a and 0. Equating the current expressions given in Eqs 3.14, 3.19, and 3.23, we find the following phase boundaries between the different current regimes ._, ^ . * 1 T a(l - a) a (Entry) a < a = — = — - JL = —— oL = Vd+l l + {d-l)a L l + (d-l)a a < f3 ( E x i t ) " < * = V l T T J ^ TTW^W °«= TW1^ <3-24> P < a 1 1 1 (Maximal) ot,B>—=—- J m a x = —T=—— omax = Vd+l m a x ( v ^ + l ) 2 tMX~ 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 wi l l see that the current and boundary densities in Eq . 3.24 are significantly more accurate those predicted by the simplest mean field theory (Eq. 3.8). In fact, the currents and boundary densities of Eq . 3.24 match our M C simulation results to within simulation error. The qualitative C h a p t e r 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 P a r t i c l e s 23 difference between the steady-state boundary densities predicted by E q . 3.8 and Eq . 3.24 arise directly from the exclusion of the individual particles. Note that while the boundary densities are accurately predicted with 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]. Al though the results in Eq . 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 in the transition between the maximal current and boundary limited regimes, while the transition between the entry and exit limited regimes exhibits a discontinuity in the first derivative of the current, and the boundary-limited to maximal current transitions have the following ^-dependent discon-tinuity in the second derivative The d dependence of Eq . 3.25 is different from that predicted using the simple mean field approach of Eq . 3.9 and is shown below to compare favourably wi th results from M C simulations. 3.3 Monte Carlo Simulations Extensive Monte Carlo simulations were performed to validate the various analytical models presented in the previous two sections. As 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]. The 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. The minimum separation restriction (3.25) C h a p t e r 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 P a r t i c l e s 24 was maintained in the simulation by setting the rate of any particle movement that yielded an inter-particle separation of less than d sites to zero. Excluding these modifications, the simulation proceeded exactly as described in Chapter 2 and Appendix A . Also, note that by modifying the B K L algorithm in 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 with larger particles, the length of the lattice was scaled wi th d to ensure that at least one thousand particles could occupy the lattice simultaneously (N = lOOOd). The simulations were run for 4 x 10 9 steps, which was 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 in [3], was used to generate the pseudorandom numbers in the simulation. The period of this generator is approximately 2 x 10 1 8 and thus is much larger than the number of simulation events. 3.3.1 Currents Figure 3.4 plots the steady-state current as a function of particle size d in all three current regimes. Whi le 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 within simulation error ( « 1%). The currents for various values of d as functions of a and (3 are shown in F ig . 3.5. From 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 in F i g . 3.5 are a combination of the results from two different equations. Equation 3.14 was used for entrance and exit rates greater than the values of a*, (3* in Eq . 3.24. For smaller values of a, (3, Eqs. 3.19 and 3.23 were used. The simulated and refined mean field curves (Eq. 3.24) are symmetric wi th 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 C h a p t e r 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 P a r t i c l e s 25 0.1 0.075 t 0.05 0.025 0.25 0.2 10.15 0.1 0.05 0.08 0.06 0.04 0.02 0 1 1 I I I I I I • M C Results — Refined Mean Field - - Simple Mean Field (a) Entry Limited Current 1 1 I I I I i i - - - _ * * ' i i i i ~ 7 1 1 - V I I I I i i i i i i i (b) Maximal Current \ '— -- _ ^ * • — ~* • -• • m i i I I I I 1 1 \ - \ . \ \ \ \ \ i ~» i_ I I I I . . I I I I i i i i i i i (c) Exit Limited Current -i i i i i i i 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 d Figure 3.4: The currents wi thin 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) The entry limited current JL as a function of d for (3 = 10 and a = 0.1. (b) The maximal current J m a x for a = (3 = 10. (c). The exit l imited current JR for a = 10 and (3 = 0.1. 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) The entry limited to maximal current profile with (3 = 10. (b) The exit l imited to maximal current profile wi th a = 10. The 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). C h a p t e r 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 P a r t i c l e s 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. The points in these plots each represent the mean of 10 M C runs of 4 x 10 9 events. The errors on the M C currents 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 M C results was within simulation error. The 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 in Eqs. 3.8 and 3.24 with those derived from M C simulations. The 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 omaxd), and the right boundary density aR are shown as functions of particle size d. The densities are the mean of 10 M C runs and have an error of approximately 3 percent. The densities or,d and crma,xd, and oR 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 oR can be readily and ac-curately computed in 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 in the entrance or exit l imited regimes are found by equating Eq . 3.13 to JL or JR in Eq . 3.24 and solving for the appropriate root o. The boundary and interior densities for all regimes are C h a p t e r 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 P a r t i c l e s 28 • oLd (a = Q.l, (3= l) t 1 1 1 1 1 1 1 1 j 1 2 3 4 5 6 7 8 9 10 d Figure 3.6: Boundary densities in each of the three current regimes. The solid lines show the predicted densities from the refined mean field approximations (Eq. 3.24). (Entry) OL a Cbulk (Exit) l + (d- l)a 1 0(1-0) d ad(l + (d - 1)0) (2d) -1 l + (d- 1)JL-y/(l + (d-l)JL)2-AdJL (2d)~ \ + (d-l)JR+ ^(i + (d-i)jRy-4djR O-R a(l - a)/0 l + (d- l)a l + (d-l)0 (Maximal) a (Vd+1) d(^/d+l)2 a Vd(Vd+l) 0(Vd + l)2' (3.26) where err, in exit and maximal current regimes are the averages of the slowly varying densities over the first d lattice sites. Recall that in 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 wi th particle size d = 3. Notice that a weak anti-correlation arises at the right boundary even in the entrance rate limited Chapter 3. Totally Asymmetric Exclusion Processes with Finite Size Particles 29 phase (inset of F ig . 3.7(a)). Similarly, there is a very small periodic behaviour in the density near the entrance in the exit rate limited regime (first inset of F ig . 3.7(c)). A more striking anti-correlation is evident near the right boundary of the exit limited lattice (second inset in F ig . 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. This 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. When 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 The 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 al l particles can only move to the right. As particles are located by their left edges, and the observed correlations are due to close-packing of the particles near the lattice boundaries, in 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 wi th finite sized particles, the bound-aries between the exit l imited, entry limited, and maximal 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 l imited phase. The following fifth order centred difference formula was then used to compute the first through third derivatives of the current data dJ dx (J(x0 - 2 A x ) - 8 ( J ( x 0 - A x ) - J(x0 + Ax)) - J(xQ + 2 A x ) ) . (3.27) Discontinuities in 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 maximal current transitions were marked by C h a p t e r 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 P a r t i c l e s 30 0.5 0.4 0.3 0.2 0.8 0.6 b~ 0.4 0.2 2 1.5 1 0.5 0 0 (a) Entry Limited (a = 0.1, (3 = 1) 0.26 0.24 0 .2^ 80 3000 (b) Maximal Current (a = 1, (3 = 1) 1 0.9 0.8 (c) Exit Limited (a = 1, (3 = 0.1) 0 10 20 2970 3000 500 1000 1500 i 2000 2500 3000 Figure 3.7: Representative scaled density (aid) profiles in the three current regimes, (a) The entry rate limited regime (a = 0.1,/? = 1). (b) The maximal current regime (a = B = 1). (c) The exit rate l imited regime (a — l,@ = 0.1). The boundary regions are shown in the insets. The strong split t ing of the densities into high density and low density branches, is a product of the close-packing of the particles in 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: The phase diagram for the T A S E P wi th d = 3 particles derived from the refined mean field results (Eq. 3.24) and the simple mean field results (Eq. 3.8). The points are the transitions estimated from Monte Carlo simulations. C h a p t e r 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 P a r t i c l e s 32 (a) (b) 0 . 8 * S 0 . 4 0 . 2 M C Results '~~ Simple Mean Field Refined Mean Field a = p (Refined Mean Field) a* (Simple Mean Field) (3* (Simple Mean Field) 5 6 d 1 0 Figure 3.9: (a) The discontinuity in the second derivative of the steady-state current with respect to a or 8 at the transition between the maximal current and the entry limited regimes. The simple mean field approach (Eq. 3.9) predicts an increasing discontinuity with d, while the refined mean field approach (Eq. 3.25) predicts a jump that decreases wi th d. (b) The critical values a* and 8* at which the boundary limited to maximal current transitions occur (Eqs. 3.8 and 3.24). peaks in the thi rd derivative. The transition point predicted by Monte Carlo was assumed to be at the value of a or 8 associated with the maximum 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 ig . 3.8. Addit ionally, and as can be seen in Fig . 3.9, the M C results for the jump in the second derivatives d 2 J / d a 2 and d 2 J / d f 3 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. Totally Asymmetric Exclusion Processes with Finite Size Particles 33 3.4 Summary We have presented mean field theories which accurately model the currents and densities of a totally asymmetric exclusion process with particles of arbitrary size. Specifically, steady-state currents and densities were computed, using both analytic, mean field ap-proaches and Monte Carlo simulations. A l l calculations show that the steady-state cur-rents qualitatively retain the three phase structure (low density, high density, maximum current) familiar from the standard (d = 1) T A S E P . Our best mean field model utilizes the canonical parti t ion 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) wi th 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. Both analytic approaches predict first order transitions in the current between the entry and exit l imited phases, and second order transitions in the current between the maximal current and the boundary l imited phases. However only the refined mean field approaches produce estimates for the discontinuities in the current derivatives that match the simulation results. The 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 in the theme of the preceding chapter, we now consider other modifications to the standard totally asymmetric exclusion process. A s described in 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 in 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 with selectivity filters and solvation zones. To improve the uti l i ty of the T A S E P in 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 with multiple movement rates focused on systems with 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 in 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 in the interior of a lattice has been investigated [17]. When the hopping rate associated with the defect site was faster than the majority of sites in the lattice there was no change J The work in this chapter is based on "Steady-state properties of a totally asymmet-ric exclusion process with periodic structure", Greg Lakatos, Tom Chou, and Anatoly Kolomeisky, Phys. Rev. E 71 (2005) 011103 C h a p t e r 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 w i t h P e r i o d i c S t r u c t u r e 35 7" T T ? W t V I Figure 4.1: The periodic dual-rate totally asymmetric exclusion model. A t every T lattice sites, the particle movement rate is p 2 . 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 T A S E P . 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 maximal 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 Model and Methods The generalization of the T A S E P investigated in this chapter includes two internal hop-ping rates, p \ and p 2 . Numbering the lattice sites starting from 1 at the far left, we assume the p 2 sites are arranged periodically with a period T and that the lattice consists of a whole number of ( p i , p i , . •. ,p2) segments (Fig. 4.1). 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 l imited regimes. The accuracy of our mean field methods are determined through comparison with results from Monte Carlo simulations. Chapter 4. A Totally Asymmetric Exclusion Process with Periodic Structure 36 4.3 Mean Field Theories 4.3.1 Simple Mean Field Methods We begin by considering the dual-rate T A S E P in the l imit of large a and B, where the behaviour of the system wi l l be determined entirely by the internal movement rates p\ and p2. Ensuring continuity of the current in a lattice wi th periodicity T we find J = p2oT(l - < T i ) = p i < 7 i ( l - a2) = . . • = p i c r T _ i ( l - oT) (4.1) In writ ing E q . 4.1 we have assumed that the densities (oi) in the lattice have the same period T as the hopping rates in the lattice. Solving E q . 4.1 for J in terms of one of the densities o"j, then maximizing J, yields an approximation to the maximal current and densities in terms of p\ and p2. Unfortunately, while it is easy to write Eq . 4.1 for arbitrary T , the solution of Eq . 4.1 yields increasingly unwieldy expressions for the current and the densities as T increases. As a result we only show the expressions for the T = 2 case: j = V\V2 2 ' * i = ^ T / - ' (4-2) /Pi a2 = ^Pl + VP2 These relations show the expected invariance under the interchange of pi &ndp2- 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 in 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 p2 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, C h a p t e r 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 37 (4.4) J = a(l- <TI) = pi<ri(l - a2) = •••= P20T{1 ~ < 7 i ) . (4.3) From Eq . 4.3, we find the following densities and currents for the T = 2 case = P i { p 2 - a)a P2(Pi + o)-api _ ap2 p2{pi + a) - api oa>2 = Oil pi. Similarly in the exit l imited case we have J = P20"/v(l ~ °~N-T+l) = PlO~N-T+l(l - C/v-T+2) = . PlCTjV-l(l - ON) = PO-N-Giving, for T = 2 Pl(p2-(3)i3 (4.5) P2(pi +/?) - / 5 p i ' ^ , 1 = 1 - - , (4-6) P2 = Pl(P2~ P) P2(Pl+P)~PlP' To determine the transitions between the maximal current and entry limited regions we equate the maximal current solutions generated by E q . 4.1 to the expressions for Ja and J/3 in Eq . 4.4 and Eq . 4.6 and solve for the transition values a* and P*. For the T = 2 case this yields a* P2\JP\ ,A a =P = r — — ( 4 . 7 ) VPI + VP^ Final ly equating the current expressions in Eq . 4.4 and E q . 4.6 we find that the transition between the entry and exit l imited regions occurs when a = 0. C h a p t e r 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 w i t h P e r i o d i c S t r u c t u r e 38 4.3.2 Refined Mean Field Methods As we wi l l see in the next section, the results of the simple mean field approach do not provide a particularly good match to the results of Monte Carlo ( 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. The 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 wi thin a single period, and neglect correlations between adjacent periods. This 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 wi l l be able to produce concise and accurate expressions for the currents and densities in all three phases of the dual-rate T A S E P . We first consider the maximal 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 ) . The pair probability P ( x j , x i + 1 ) is the probability of finding a p i site with occupancy X ; , followed by a p 2 site with occupancy xi+i. The time evolution of the occupancy state of any two adjacent sites in a T A S E P wi l l depend on the two sites themselves along with 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 Kdt' = - p 2 ( P ( 0 , l , 0 , 0 ) + P ( l , l , 0 , 0 ) ) + p 2 ( P ( 0 , l , 0 , 0 ) + P ( 0 , l , 0 , l ) ) . (4.8) Now assume that each pair of (p i ,p 2 ) sites behaves like a statistically independent unit, uncorrelated with the other (p\,p2) periods in the lattice. Using these assumptions the probabilities in E q . 4.8 can be decomposed into products of pair probabilities, P ( X J , xi+i, xi+2, xi+3) « P ( x i , x i + i ) P ( x i + 2 , xi+3). This yields the following equation for P ( 0 , 0) in the steady-state: Chapter 4. A Totally Asymmetric Exclusion Process with Periodic Structure 39 = p 2 ( P ( 0 , l ) 2 - P ( 0 , 0 ) P ( l , l ) ) = 0. ( 4 9 ) Additionally, given the definition of P(xi,Xi+i), the following relations hold 01 = P(1,0) + P(1,1), 02 = P(0,1) + P(1,1), (4.10) P(0,0) = 1 - P ( 0 , 1 ) - P ( 1 , 0 ) - P ( 1 , 1 ) . Using Eq. 4.10 in Eq. 4.9, produces the relation 7 ^ ' . (4.11) 1 + (cr 2 - O-i) An equation for (5(1,0), the probability of having an occupied p2 site followed by an unoccupied pi site, can be found from Eq. 4.11 by interchanging p1 and p 2 , as well as o\ and o2. Applying the current continuity condition PiP(l ,0) = p2Q(l,0), results in the relation o-i(1- cry) - (02 ~ oif _ o2{l - Qi) - ( " i - 02? 1 + ( fJ 2 - <7l) 1 + (<Ti - <72) Now consider the substitution o2 — 0\ = /c and assume pi to be the larger of the two rates. Thus oi < o2 and k > 0. Solving Eq. 4.12 for o\ yields ol = (2k(Pl + p2) - 2{Pl - p2)Yl ( p 2 ( l - k2) - Pl(k - l ) 2 ) + V ( l - k2){k(Pl +p2) - (Pi - P2)){Sk(Pl + p2) - {Pl-p2)) (4-13) Equation 4.13 shows that o\ is real when — 1 < k < 3 ^ ^ ) or < k < 1. Substituting Eq. 4.13 into Eq. 4.12 gives J = -, kflVJ , >, which shows J < 0 for k > J > 0 1 1 fc> ( P l - P 2 ) - « ( P 1 + P 2 ) ' P 1 + P 2 ' for 0 < k < 2i-jE2, and § = P I P ^ - P ' ) > 0 V k. As J must be positive, and P 1 + P 2 ' dk ( f c ( p 2 + p i ) + ( P 2 - p i ) 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 0"2 J = Pi + 2p2 3(pi + p 2 ) ' 2p i + P2 3(pi+P2)' ( 4- 1 4) P1P2 2 ( p i + p 2 ) ' Note that the solution is invariant under the exchange ( c r ^ p i ) for ( c r 2 , p 2 ) , and reproduces the standard TASEP results in the limit pi = p2: 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 l - ^ ^ l = n ^ ^ - v (4-15) 1 S(Pl+p2)J 2(Pl+p2) v ; 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 = 3PlP2 (4 17) P 2 (2pi+p 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 Exclusion Process with Periodic Structure 41 rfP(0,0) dt riP(0,l) Jt dP(l,0) dt dP(l, 1) = -aP(0 , 0) + p 2P(0,1)(P(0,0) + P(0,1)), = -c*P(0,1) -p 2P(0,1)(P(0,0) + P(0,1)) +PiP(l,0), = aP(0,0) + p 2 P ( l , 1)(P(0,0) + P(0,1)) -P iP( l ,0 ) , = - p 2 P ( l , 1)(P(0,0) + P(0,1)) + « P ( 0 , 1 ) . (4.18) dt Similarly, in the exit limited phase we find the following master equation for the last two sites in the lattice rfP(0,0) dt dP(0,l) dt rfP(l,0) dt dP( l , l ) /5P(0, l)-p 2 P(0,0)(P(0, l) + P( l , l ) ) , -/3P(0,1) - p 2P(0,1)(P(0,1) + P ( l , 1)) +PiP(l,0), ( 4 . 1 9 ) /?P(1,1)+ P 2P(0,0)(P(0,1) + P(1,1)) -P iP( l ,0 ) , - -/3P(1,1) +p 2P(0,1)(P(0,1) + P ( l , 1)). dt Equations 4.18 and 4.19 assume that the occupancy probabilities P(xi,xi+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 TASEP 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 TASEP (c. f. Fig. 4.6). Applying the normalization condition on the probabilities P(xi,xi+i) and solving Eqs. 4.18 in the steady-state gives two sets of solutions for the C h a p t e r 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 w i t h P e r i o d i c S t r u c t u r e 42 current and bulk densities in the entry limited phase. We choose the solution which gives o~\ —• 0, c 2 —> 0 in the limit a —> 0. This yields the following results for the current and bulk densities in the entry limited phase « \P1P2 - a ( p i + a) + y/4pla(p2 - a) + (pi(p2 - <*) + a 2 ) 2 ) 2p2(Pi + a) _ ap 2 + (a + pi)(a + p2) - y / b p j a i p i - a) + (pi(P2 - a) + a 2 ) 2 (4.20) a ' X 2p2(p! + a) a C a , 2 = — • P2 Similarly, solving Eqs. 4.19 and choosing the solution that gives u\ —> 1, cr2 —> 1 in the limit P —+ 0, produces the following approximations to the current and bulk densities in the exit-limited phase P [PlP2 - P(Pl +P)+ V4P(P(P2 ~P) + (Pl(P2 ~P)+ P2)2 J P = 2p2(PL+P) ' *M = 1 - A ( 4- 2 1) P2 PXP2 ~ P(Pl +P) + VMPJP2 ~P) + ( P l ( P 2 ~ P ) + W ^ 2p 2( P l+/3) 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 — cra>1) = o p ^ , and (1 — <7 a ) 2) = o/3ti 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 TASEP 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 C h a p t e r 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 w i t h P e r i o d i c S t r u c t u r e 43 4 _ °eff P2 •'2 P, 1=1 1=2 r=3 1 = 4 L = 4 Segment eff 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 /3eff = p 2 ( l — ° i ) o n the right. 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 (pi ,pi , . . . , p2) periods (Fig 4.2). Assuming the bulk densities in the TASEP lattice are periodic and inherit the same periodicity as the movement rates, we define an effective injection rate (cteff) and an effective extraction rate (/?e~) for this finite segment. Specifically, we set "eff = P2OL and /3eff = p 2 ( l — (7i). Using these rates and the known values of p\ and ~>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 M P ( t ) . (4.22) Here M is the matrix of transition rates that describe the TASEP and P(t) is a vector of the time dependent probabilities of finding the TASEP lattice in a particular occupation state. The transition matrix M is typically sparse, and even for a small segment it is also large (2L x 2L). To construct M efficiently it is desirable to avoid having to test the 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 t h entry of f equal to the hopping rate of the i t h lattice site in the finite segment. The hopping rate of the last site in the finite segment however is set to /?efr> so that rL — B e f i . With the finite 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: The steps in 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 wi th the corresponding decimal value (i.e. lattice occupancy O i l —> state 3). The 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 in 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 in the figure). Cal l ing 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 ina l ly 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 matrix is (see F i g . 4.3) 1. L a b e l t he o c c u p a n c y s tates o f the finite segment Label each occupancy state of the finite-segment wi th a number determined by treating the occupancy of the lattice as an integer value expressed in base-2. For example, an L = 3 lattice with only the last site filled (i.e. an occupancy of {001}), is said to be in state 1, while the same lattice with the last two sites filled (an occupancy of { O i l } ) is said to be in state 3. 2. Sepa ra t e t h e s ta tes i n t o t w o g roups Separate the states labelled in 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. Enumerate 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 2 ; the rate associated with the second 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 aerf. For example, if we are working 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 aeff. 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 TASEP lattice. This iterative procedure proceeds as follows 1. Bui ld the Master Equation Given an initial guess for the values of (<TI, CTL), use the recursive algorithm to pro-duce the transition matrix (M(<7i, OLIVIIVI)) f ° r the finite-segment master equation. The finite segment densities (OI,OL) a r e introduced into M solely through the ef-fective injection rate ctes = P2&L, and the effective extraction rate /?eff = p 2 ( l ~~ 2. Solve the Master Equation To find the steady-state currents and densities, compute the eigenvector (V(ci,o2,Pi)P2 of M(CTI , cr i ,p 1 ,p 2 ) associated with the zero eigenvalue M(o1,o2,Pl,p2)V = 0. (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. Compute 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 Structure 47 E v> (4.25) 4. Start a new iteration or finish Compare the new density estimates (ol,o*L), to the previous estimates (OI,OL)- If |cr* — o-!| < e and \o*L — OL\ < e, then the iteration has reached a fixed point and the procedure ends. Here e is an arbitrary convergence parameter (set to 10~4 in this study). If the convergence criteria are not met, return step 1, setting o\ = cr*, and OL = o*L. Note that as only the effective injection and extraction rates change, 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 For convenience we refer to the use of these two algorithms to compute the steady-state properties of a TASEP 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 Qferf = a or (3e$ = (3 as appropriate, and applying the density self-consistency condition at the end of the finite-segment that lies in the interior of the TASEP 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 exponen-tial 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 J = a e f f ( l - o-i) = P2<yL{l - oi). (4.26) 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 TASEP results (Pl = p2) by less than half a percent. As a result, unless otherwise noted, we used lattices containing 1000 periods for all our simulations. The simulations were run for 7 x 108 steps, which was sufficient to reproduce the known TASEP 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 Exclusion Process with Periodic Structure 49 (a) 0.25 0.2 ^ 0.15 0.1 0.05 0 (b) 0.25 0.2 ^ 0.15 0.1 0.05 Q • M C Results — Simple mean field --• Refined Mean Field — L = 12 FSMF 0.4 0.6 " 2 Figure 4.4: (a) The 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 maximal current phase (a = 10, (5 = 10). (b) Max ima l currents for dual-rate T A S E P s of various periods. The F S M F estimates shown by solid lines were produced using a single period of the lattice as the finite segment. Whi le 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 with 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 p2/Pi ~ 1. In this limit the model is essentially the normal single rate TASEP, and the expected pair correlations between adjacent p\ and p2 sites are small. 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 p2. The F S M F method clearly performs better for smaller values of T. For larger T, the method underestimates the Monte Carlo currents for 0.1 < p2 < 1.0 and overestimates for very small values of p2. 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 mean-field coupling of the finite segment to the remainder of the lattice. The fact that the finite segment underestimates the Monte Carlo current for p2 > 0.1 suggests the existence of an anti-correlation between the p2 site of one period and the first pi site of the next. Additionally, the fact that the FSMF 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 TASEP 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 Exclusion Process with Periodic Structure 51 0 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 TASEP 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 indistin-guishable, (b) Current profiles for dual-rate TASEPs of various periods with p 2 = 0.25. The F S M F results (lines) were produced using an L = T FSMF approach (/? = 10). Chapter 4. A Totally Asymmetric Exclusion Process with Periodic Structure 52 0.6 W) 0.5 0.4 0.3 0.2 0.1 0.8 .0.6 0.4 0.2 0.9 0.8 0.7^ 0.6 \ 0.5 P 0.4 Entry Limited (a - 0.05, (3 = 10) (b) Maximal Current (a = 10, (3 = 10) [(c) Exit Limited (a = 10, (3 = 0.05) 0 500 1000 i 1500 2000 Figure 4.6: Monte Carlo density profiles for a T — 2 TASEP with p2 = 0.1. Comparison with Fig. 2.1 shows that the profiles in the T = 2 TASEP retain the same qualitative character as the profiles in the standard TASEP. 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. 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 Pi and p2 sites in the centre of the TASEP lattice. The mean density at the two classes of sites in the TASEP 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 Exclusion Process with Periodic Structure 53 ef 1 0.8 ,0.6 0.4 0.2 Q i i i i i i —"T —1 — 1—I — 1 I I I I Refined Max Current MF - - Simple MF - - L = 4 FSMF / P / / i _l 1 I I I I ..I. .1 . 1. i. I_t-0 0.2 0.4 0.6 Pi 0.8 Figure 4.7: Maximal current phase (a = 10,/? = 10) densities for a T = 2 TASEP 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 p2 where the density correlation between adjacent sites is expected to be large. periodicity in all three phases. While the FSMF 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 p2 = 0.2. In the maximal current phase the profiles in any one period are approximately linear which is similar to the behaviour of a TASEP 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 p2 ~ 0.5. The exit-limited results are related 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 Exclusion Process with Periodic Structure 54 -10 -5 0 5 10 i Figure 4.8: The bulk density profiles for the T = 5 and T = 9 dual-rate TASEPs with p2 = 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 FSMF (lines) results. The M C results are the central 5 periods (T = 5), and central 3 periods (T = 9) of the lattice. 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 TASEP 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 TASEP lattice. The F S M F results are indistinguishable from those of Eq. 4.20. Although the solution for o2 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 o2 profiles predicted by the two mean field methods at small values of p2. Chapter 4. A Totally Asymmetric Exclusion Process with Periodic Structure 56 4.4.3 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 a2 in the centre of the lattice, and locating any clear discontinuities in the derivatives. The dual-rate TASEP retains the general form of the standard TASEP 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 (p2 in our example). Physically, 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 p2 is then expected, as decreasing p2 slows the motion of the particles in the interior of the lattice. Despite the varying degrees of success in predicting accurate steady state currents and densities, all three mean field approaches predict the phase boundaries within approxi-mately 20% error. Note that the phase transitions predicted by the boundary-limited re-sults 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 p2 = 0.1 is held constant. The major effect of increasing T at a fixed p2 is to increase the values of a* = ft*. 4.5 Other Periodic Arrangements The dual-rate TASEP model investigated in the previous sections was built on lattices composed of an integer number of {p\,Pi,Pi, • • • ,p2} periods. While a useful extension to 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 C h a p t e r 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 w i t h P e r i o d i c S t r u c t u r e 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 two-rate TASEP at various periodicities, with p 2 = 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. Chapter 4. A Totally Asymmetric Exclusion Process with Periodic Structure 5 8 (a; 6 = 1 a (b) 8 = 2 a ( c ) 5 = 3 a 8 9 '1 Pi -T—-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 TASEP 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 TASEP model however, the issue of particle-hole symmetry arises. To investigate the conditions under which a periodic dual-rate TASEP would preserve particle-hole symmetry, first define a vector f consisting of all the movement rates in the TASEP 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 TASEP lattice. Now consider the probability P(a, 6, r, xi, x2) • • •, XN) of finding the TASEP lattice in an occupancy state xi,x2, • • • ,x.jv with injection rate a, extraction rate /?, and hopping rates r. Here x.h = 1 if lattice site i is occupied, Xi = 0 otherwise. Particle-hole symmetry then requires P(a,f3,f,{xi,x2,...,xN}) = P(B,a,f,{xN,xN-i,...,xl}) (4.27) where = 0 if Xj = 1, and xx — 1 if xx = 0. For Eq. 4.27 to hold, the rate of the particle movement that takes the TASEP lattice from a occupancy state s\ = {xi ,x 2 , . . . to an occupancy state s2 must equal the rate of the particle movement that takes the lattice from state S\ = {xN,xN_i,... ,x{\ to occupancy state s2. This places constraints on the 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 2 - By definition this movement would occur at rate r^. In lattice state si the equivalent 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, must equal r u - h - Now consider the case where = p 2 . Then by definition of the dual rate 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 TASEP 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 2 , consider the case when k = m T + 7 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 = r(u!-m)T+6+(5-y) (4.29) 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 2 under the symmetry exchange operation. An interesting consequence of the restriction on is that a lattice with the first p 2 site located at site T and containing a whole number of periods preserves symmetry regardless of the value of T, while a lattice with the first p 2 site at lattice site 1 and containing a whole number of periods violates symmetry if T ^ 2. Figure 4.12 shows the phase diagrams for various symmetric and non-symmetric re-alizations of the dual-rate TASEP with T = 3 and p 2 = 0.25. The results displayed in C h a p t e r 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 w i t h P e r i o d i c S t r u c t u r e 60 ( a ) 1 0 . 8 0 . 6 CO. 0 . 4 0 . 2 ft f 4 - 5 = 1,N = 1 2 0 5 - 5 = 2,N = 1 2 0 4 - 5 = 3,N = 1 2 0 3 • 5 = \,N = 1 2 0 3 • 5 = 2,N = 1 2 0 3 • 5 = 3,N = 1 2 0 4 0 . 2 0 . 4 cx 0 . 6 0 . 8 Figure 4.12: Phase diagrams for the T = 3 dual-rate TASEP for various values of 8 and N derived from Monte Carlo simulations (p2 — 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) TASEP and the (5 = 2,N = 1204) TASEP 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*. 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 TASEP 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 TASEP 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 TASEP and the symmetric version of the T = 3 dual-rate TASEP with (5 = 1,N = 1205) both place the first defect at lattice site 1. The two realizations of the TASEP 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 TASEP is greater than the value of p2 [a* « 0.27, p2 = 0.25). When we consider that the first lattice site of the TASEP is a defect site, then the effective injection rate into the remainder of the TASEP lattice cannot exceed the defect rate. As a result the TASEP will be in an entry limited phase until a at least matches p2. 4.6 Summary We have developed approximate, mean field, methods to compute the steady-state cur-rent and densities of totally asymmetric exclusion processes involving two internal hopping rates. Additionally, we have simulated the two-rate TASEP and have explored its phase diagram. We find that the dual-rate TASEP retains the three phases found in the stan-dard TASEP 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 ap-proach 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 al-ternative 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 TASEP possesses particle-hole symmetry when the number of sites in the TASEP lattice equals mT + 25 for an integer m > 0. When the size condition was not met, the phase diagram of the dual rate TASEP 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 TASEP lattice, while the value of B* is independently determined by the position of the last defect. 63 Chapter 5 Clustered Bottlenecks in mRNA Translation and Protein Synthesis 1 5.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 TASEP to the topic of biological protein synthesis. During protein synthesis, ribosome molecules bind at one end (the 5' end) of a messenger RNA (mRNA), translate along the mRNA 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 RNA (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 RNA must bind to the ribosome mRNA 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), ATA (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 tRNA molecules complementary to these defect codons can be as much as 30 times lower than the concentration of the more abundant tRNA species [8]. x The 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 C h a p t e r 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 unidirec-tionally along the mRNA as tRNAs (not shown) deliver the appropriate amino acid to the growing protein. Codons with low concentrations of corresponding tRNA ("slow-codons") result in bottlenecks that locally slow ribosome mo-tion. 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 p2 < pi, 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]. Statistics indicate a higher occurrence of rare codons near the 5' binding site of E. coli 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 tRNA 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 DNA 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 densi-ties 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-C h a p t e r 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 6 6 (a) Pi 7> • L = 5 ^ k= 1 (b) 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 sepa-rated 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 TASEP onto the mRNA translation process, the TASEP lattice repre-sents the mRNA, and each TASEP particle represents a translocating ribosome. Each time a TASEP particle makes a hop, a codon is read, and a tRNA delivers its amino acid to the growing polypeptide. When a TASEP 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 TASEP lattices approximately the same size as those investigated in the pre-ceding two chapters. The specific variant of the TASEP 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 TASEP 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 TASEP will behave like a standard defect free TASEP, we can define CT_ to be the density in the region of uniform Chapter 5. Clustered Bottlenecks in mRNA Translation and Protein Synthesis 67 hopping rates to the left of the defect cluster. Similarly we define o+ 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 TASEP 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 Pi in the uniform regions has been normalized to one. This implies that either cr_ = o+ or cr_ = (1 — o~+). Given that the T A S E P lattice contains defects, we 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 aeff = (3es = (1 — o + ). 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 TASEP lattice should be an extremely good approximation. With these assumptions, the FSMF 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 FSMF 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 TASEP behaves like a standard defect-free TASEP, 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 tRNA 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. C h a p t e r 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 6 8 0.3 0.25 ^ 0.2 ^ 0 .15 0.1 0.05 0 i i i i | i i i i | i 1 1 1 | 1 I I 1 | 1 T [ [ (a) Current, single defect / , , , , I , , , , I , , , , I , , • - • M C Simulation - - L = 1 F S M F L = 3 F S M F - L = 5 F S M F 0 0.2 0.4 0.6 0.8 1 cTO.75 3 0.25 ,1 1 1 1 1 1 1 1 1 1 1 1 1 1 ; \ (b) Current, m contiguous defects i = 0.01 — P 2 = = 0 .10 • - * p 2 = = 0 .20 ~ P 2 = = 0 .30 1•—. : ^^^^ _ - A P 2 = = 0.50 ~ *• - -' 1 " t l t 1 l 1 1 l l l 1 l l i i 1 2 3 4 5 * 0 . 5 , - _ _ _ _ J ^ m Figure 5.3: (a) Comparison of steady-state currents Ji(p 2) for a single defect derived from M C simulations with those obtained from a //-segment FSMF 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 ~>2 is changed. Thus for clarity, the currents in the figure have been normalized by the slow-codon hopping rate (p2). Also note that for small ~>2, J 2 / p 2 ^ 1/2 consistent with Eq. 5.3. Chapter 5. Clustered Bottlenecks in mRNA Translation and Protein Synthesis 69 5.3 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 finite-segment mean field approximation to the ribosome current to lowest order in the defect codon translation rate p2. 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 aeff = (Pi/p 2)(l — cr+) and the effective extraction rate to /?eff = (1 — c + ) . Making use of the exact matrix product solution from [15] for a TASEP of length m with an internal hopping rate of 1, the finite segment expression for J is j, _ Rm-i(P~j) - Rm-i(u~i) , . J ' P 2 ~ D (ro-l\ c / - I N • \PA) Where Rm is given in Eq. 2.2. Now consider the terms in J that depend on l/ae ff m + l = r ( 1 _ } -„ t—1 m!(m — n + 1) r= (5.2) m + l . . . . - v I E (n — l)(2m — n)\, . . „ . m , M l l J - , + l ) ! f e / a l " - > + ° ' ° ' » ' n=2 v ' As o+ is the density downstream from the block of contiguous defects, in the limit of small p2, o+ must approach 0. Thus the lowest order term in the power series expansion of <7+ must be at least linear in p2. The sum in Eq. 5.2 starts from 2 and thus Rm(l/aeg) and Z?m_i(l/a efj) are both of at least order 2 in p2 and can be ignored in the computation of the lowest order term in J. Then, lowest order in p2 Chapter 5. Clustered Bottlenecks in mRNA Translation and Protein Synthesis 70 E ( n - l ) ( 2 m - (n + 2)! i ^ - ! ^ 1 ) ^ ( m - l ) ! ( m - n ^ U & f f ) ^ ( n - l ) ( 2 m - r c ) ! fm)!(m + 1 — ri)\ (5.3) n=2 m + l As expected in the limit p2 —* 0, J approaches 0. Similarly, in the limit of large m, the current predicted by Eq. 5.3 approaches J = P 2 / 4 , the current of a long uniform lattice with hopping rates p2 in the maximal current regime. Figure 5.3(b) plots Jm(p2)/p2 for various p2 as a function of m defects using a large L = 14 finite segment. Note that as the currents for a wide range of p2 values are displayed, all currents in Fig. 5.3(b) are normalized by the slow-codon hopping rate p2. For all p2 < 1, Jm(P2)/P2 approaches 1/4 as m —> N. For strong defect codons (p2 < 0.3), the largest decrease in Jm{p2) occurs as m goes 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 J2(p2,k) across a lattice containing two defects spaced k sites apart. Results for k < 5 were computed using an L — 14 FSMF 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 FSMF and M C simulations in Fig. 5.4(b) for a pair of defects (p2 = 0.15). Also note that Fig. 5.4(b) shows that the density profile between the two defects is approximately linear, similar to the known behaviour of the standard TASEP 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 aeff = ArfF = P 2 ( l ~ "+)• Applying the exact matrix-product method solution from [15] to the finite segment gives J = lim - 1 M (5.4) Rk(8-^) - Rk(a-J) As the limit in Eq. 5.4 is an indeterminate form, we apply l'Hopital's rule to find ™ 2 . (5.5) /c+i &=l/<*ef  y^7n(fc)^eff n=2 (fc+l)-n In Eq. 5.5 7n(/c) = n^~fc+i2^)^! • As o;efr ~ P 2 the lowest order term in the numerator of Eq. 5.5 must by linear in p2, while the lowest order term in the denominator must be independent of p2. So that, to lowest order in p2, j ^ jk(k - l)o!ef f _ kp2(l - o+) lk+i{k) fc + l (5.6) _ kp2 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 p2. 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. Figure 5.5(a) plots the variation in J2GD2, k)/J\(p2) as a function of k for various p2, and shows that in the limit 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 TASEP 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. C h a p t e r 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 (p2 = 0.15) of various spacings k from FSMF 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. Chapter 5. Clustered Bottlenecks in mRNA Translation and Protein Synthesis 73 1 1 1 1 1 1 0 5 10 15 20 25 30 35 40 m 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 = 300 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 2 , we can estimate the expected current for m randomly 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 ~ / »rs (N -(m - l)(k- 1)\ Zk(m,N)=l 1 JK J). 5.7 \ m J The probability that at least one pair of defects is separated by exactly k sites is then Q k ( m , N ) = { Z k - Z k + l ) . (5.8) ^ i As 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,p2)RAK < Qk{m,N)J2(p2,k). (5.9) 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. The current J(m,N = 300)RAN normalized by J\{p2) is shown in Fig. 5.5(b). For very small m, the current is approximately that of a single defect. 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 Summary We 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 DNA Hairpin Unzipping 1 6.1 Background Over the last two decades, DNA biotechnology and nanofabrication have been trans-formed 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 mem-branes [37-41], as well as investigations of force-pulling mechanisms for inducing DNA dehybridization ([48-50]). Here, dehybridization is defined as the separation of the two stands of double-stranded DNA, accomplished by breaking the bonds between complimen-tary DNA bases. While the process of single stranded DNA transport through membrane channels has been extensively studied [44-46, 49], only recently has the process of driven double stranded DNA dehybridization by transmembrane nanopores been investigated [51-53, 85]. In this chapter, the dehybridization of DNA resulting from the pulling of DNA 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 po-tential was applied to drive the translocation of DNA 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 DNA strand was inside the trans-membrane channel, the channel was blocked, and the ionic current fell below the open channel value. When the DNA strand exited the channel, open-channel ionic currents resumed and were measured. The experiments measured the distribution of first passage lThe work in this chapter is based on "First Passage times of driven D N A hairpin unzipping", Greg Lakatos, Tom Chou, Birger Bergersen, and Gren N . Patey, Phys. Biol. 2 (2005) 1-9 Chapter 6. First Passage Times of Driven DNA Hairpin Unzipping 77 times for the escape of the DNA 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 DNA hairpins with single defects (non-homologous base pairs). Since there are effectively fewer base pairs to break, DNA 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 DNA 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 DNA base. The energy barrier should be ap-proximately the difference in free energy between hybridized and fully dehybridized DNA. The DNA 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 two-dimensional multi-barrier stochastic model. The DNA hybridization energetics are mod-eled using both a simple free energy model where each G C base pair contributes 3kBT, and each AT 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 DNA 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 DNA escape. 6.2 Model We model the process of DNA extraction through a membrane pore using the zipper model displayed in Figs. 6.1 and 6.2. The model assumes that the DNA 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 DNA strand are enumerated starting from the 3' end, index n 0 indicates the position of the DNA strand relative to the start of the transmembrane region of the cv-hemolysin pore. As the DNA strand can slide in reverse, no can take on negative values. Similarly, index nx equals the number of hairpin base pairs that have been dehybridized (See Fig. 6.1(b)). Using this nomenclature, the DNA hairpin is fully dehybridized when n\ = N. As the dehybridized DNA 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 0 = —NA absorbing. Forward escape (escape into the trans chamber) can only occur after the DNA hairpin has completely dehybridized, and so we also make state (n 0 = 2N + NLOOP + NL, nx — N) absorbing. Defining P(w,v,t) to be the probability of being in state (n 0 = w,ni = v) at time t, the system then evolves according to a master equation In Eq. 6.1 the expression L(i,j\w,v) is the rate of transition from the state with (n 0 = i^m = j) to the state (n 0 = w,ni = v). We restrict the set of allowable state transitions to include only four moves: sliding of the entire DNA strand in the trans direction (no —> (n 0 + 1)), sliding in the cis direction (n 0 —> (n 0 - 1)), dehybridizing a DNA base pair (ni —» (ni +1)), or rehybridizing a DNA base pair (m —> {nx - 1)). In all cases we assume dP(w, v, t) ~ dt ^L(i,j\w,v)P(i,j,t) - P(w,v,t)^2L(w,v\i,j). (6.1) C h a p t e r 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 U n z i p p i n g 79 Parameter Description V Number of base pairs in the DNA hairpin NA Number of bases in the poly-A tail attached to the 3' end of the DNA hairpin (50). NL Number of bases that can simultaneously occupy the transmembrane part of the pore (12) Nv Number of bases that can simultaneously occupy the vestibule part of the pore (16) A^ loop Number of bases in the hairpin loop (6) MR The 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 DNA strand, and are not used as fitting parameters. Figure 6.1: (a) Division of the a-hemolysin channel into the cis, vestibule, transmem-brane, and trans regions. The DNA hairpin is inserted into the channel from the cis region. The number of DNA bases that can fit into the vestibule re-gion is set to NV = 16, while the number of DNA bases than can fit into the transmembrane region is NL = 12. (b) An example of the DNA hairpin translocation experiment. The process of DNA 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. Chapter 6. First Passage Times of Driven DNA Hairpin Unzipping 80 CD C A A A A A A A A A A A 3 ' 7 6 5 4 3 2 1 Base Index Figure 6.2: The enumeration scheme for DNA base pairs in the hairpin and the associated state space for the 2D stochastic model. In this example N = 7, NA = 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 n0 have attempt frequency p.. States with nQ — 2N + N\OOP + NL (a DNA strand pulled completely into the trans chamber), and n 0 = —NA (a DNA strand pulled completely into the cis chamber) are absorbing. Chapter 6. First Passage Times of Driven DNA Hairpin Unzipping 81 that L(i,j\w,v) has the form L(i,j\w,v) = f ( i J \ 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 to state (w,v), while AG(i, j\w,v) is the free energy difference between states (w,v) and To de-fine AG(i, j\w,v) 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 DNA 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 DNA strand, from each of the four regions. In the trans region, all the DNA 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 DNA bases in the trans region to the change in free energy upon going from state (i,j) to state (w,v) is then In Eq. 6.3, Ntrans(no,ni) is the number of DNA bases in the trans chamber in state (n 0,ni), and b is the persistence length of single stranded D N A (approximately two inter-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 where A 7 " T M P (n 0 , n\) is the number of single-stranded DNA (ssDNA) bases stored in the (6.3) A G T M P ( « , j\w,v) = kBT(NTMP(w,v) - NTMP(i,j)) (6.4) Chapter 6. First Passage Times of Driven DNA Hairpin Unzipping 82 transmembrane part of the pore in state (n 0 ,n i) . Turning to the vestibule region we 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 AGvest(iJ\w,v) = svT(Nvest(w,v) - Nvest(iJ)), (6.5) where A\/ e stibuie(^o, n \ ) 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 sv > 0 is the entropic cost for storing a DNA base in the vestibule region. While the hairpin is in the vestibule region and sv is non-zero, the free energy cost to open a base pair is increased by 2svT and the hairpin is stabilized. While we leave sv as a fitting parameter, the confinement free energy can be estimated from the entropy cost per persistence length to confine a polymer to a narrow tube: sv ~ k,B{l/D)i [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 DNA as a polymer tethered to the membrane. With this free energy, the contribution to A G from the cis region is ([42, 44, 47]) AGcis(i,j\w,v) = ^kBTln where Ncis(n0, 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 DNA hairpin. Next, we introduce the enthalpic contribution to the free energy from the change in the electrostatic potential associated with sliding a charged DNA strand Ncis(w,v) Nds(i,j) (6.6) AGel(i,j\w, v) = U(w,v, V) - U(i,j, V). (6.7) Chapter 6. First Passage Times of Driven DNA Hairpin Unzipping 83 Here U(n0, nX) V) is the potential energy of the D N A strand and is given by NA+n0 U(n0,ni,V)= £ qcP(i,V) i=lb (6.1 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 l0 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 width of the transmembrane part of the channel [44], and that the voltage profile in the transmembrane channel is linear 0 if i < 0 \ 1/2)^- ifO<i<NL V if i > NL (6.9) To complete A G we include the free energy cost associated wi th 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 and 3kBT, 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 in 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)) + AGe](i,j\w,v)+ AGcis(i,j\w,v) + AGvest(i,j\w,v) + AGTMp(i,j\w,v) + AG trans (hj\w,v). (6.10) Here g(ni) is the free energy of a hairpin in bulk solution wi th 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 with 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 DNA 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 DNA, 7 is a measure of the stiffness of the DNA 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 5s _ 1. Using the value of C given in [44] for a persistence length of single stranded DNA confined in the transmembrane region gives an estimate of p ~ 106s~"x. We constrain p to lie within the range of these two estimates. To find an estimate for u>, we use a Morse-like potential to model the hydrogen bonds between two DNA bases [48]. A calculation of the Kramer's escape rate to open a single base then gives a rough estimate of u> » 1 x 10 1 2s _ 1 , consistent with previous studies of DNA dehybridization [48]. While the DNA 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 DNA 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 DNA translo-cation and dehybridization process. For example, by assuming the DNA hairpin dehy-bridizes only from the free end, the effects of bubble formation are ignored by our model. Here, bubbles refer to dehybridized regions of DNA that are bounded by at least one Chapter 6. First Passage Times of Driven DNA Hairpin Unzipping 85 Parameter Description At tempt frequency for sliding transitions (n0 —> n0 ± 1) u Attempt frequency for hybridization transitions (rii —> n i ± 1) q Mean partial charge per D N A base Sy Entropic penalty for confining a s s D N A base to the vestibule region Table 6.2: List of fitting parameters used in the two dimensional stochastic model. hybridized base on both sides. Under physiologic conditions, structural fluctuations in 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 wi th base-pairing regions approximately 10 bases in length, and thus considerably smaller than the typical bub-ble size found in D N A at experimental conditions. Addit ionally, recent experiments [90] have studied the formation of bubbles in small, atypically bubble-prone po ly -AT D N A molecules approximately 18 bases in length. Here, the formation of D N A bubbles approx-imately 2 — 10 bases in length was observed in the D N A strand. However, Altan-Bonnet et al. [90] were able to measure the typical lifetime of these small bubbles to be ap-proximately 50^,s. This lifetime is typically much smaller than the mean time for the dehybridization and extraction of small D N A strands through a-hemolysin at moder-ate 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 hybridiza-tion region of the hairpin is offset by one or more bases (i.e., in the arrangement of F ig . 6.2, complement base 2 hybridizes with base 5). Under the assumption that the double-stranded 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 l ip id 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 DNA 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 DNA 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 DNA coil to diffuse this distance. When fully dehybridized, the strands we consider are 76 bases long and have relaxation times ~ 10_ 9s, considerably smaller than the ~ lCT 6s typically required to translocate the DNA strand by a single base. 6.3 Results and Discussion Equation 6.1 is used to compute the thermally averaged time for the DNA hairpin to first reach either set of absorbing states (n 0 = —NA,UI) or (2N + + A ^ n i = N). To do this, the (n0 = — NA, n{) and the (n0 = 2N + Afi o o p + NL,ni = N) states are rendered completely absorbing by setting P(-NA, ni)(t) = 0 and P(2N + A^ l o o p + NL, N) = 0 for all t. Equation 6.1 is then dP(w, v, t) dt ^(2N+Nloop+NL,N) = /~^L{i,j\w,v)P(i,j,t) - P(w,v,t)^2L(w,v\i,j). (6.12) Defining A(t\(n'Q,n[)) to be the probability that the system reaches an absorbing state between time t and t + dt, given that it started in state (n^n^) yields M t \ ( « ) ) = d_ "dt E yt(2N+Nloop+NL,N) = L(2N + Nloop + NL-1, N\2N + Nloop + NL,N)x P(2N + Nloop + NL-l,N,t)+ N J2 L{-{NA ~ l)M - NA,m)P(-(NA - l),nut). n i = 0 (6.13) Chapter 6. First Passage Times of Driven DNA Hairpin Unzipping 87 Name Sequence M F O L D 2-3 Defect Free G C T C T G T T G C T C T C T C G C A A C A G A G C -31.9/cBr -26.0kBT Defect-Bearing G C T C T G T T G C T C T C T C G C A A C T G A G C -27kBT -24MBT Table 6.3: DNA 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'0,n[)) computed from the solution to Eq. 6.12, the mean first passage time to an extracted state is /•oo e (n i ,n ' 2 ) = / tA(t\(n'1}n'2))dt (6.14) Jo 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 e-pAG(o,o\i,j) • (b-ib) where AG(0, 0|i, j) is the free energy of state relative to state (0,0). In the experi-ments of Mathe et a/., the DNA 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 defect-free 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 Unzipping 88 Energy Surface M s - 1 ) Q Sy M F O L D 1.4 x 106 1.5 x 101 3 0.39e 0.10kB 2-3 7.8 x 105 5.3 x 101 0 0.30e omkB Table 6.4: Parameters from fitting to the experimental results for the mean escape time of DNA 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.9kBT difference in hybridization free energy between the two strands may be an over-estimate. The value of the vestibule confinement entropy (sv) is also very small, suggesting that the DNA in the vestibule region has considerably more conformational flexibility than DNA stored in the transmembrane region of the channel. We also note that the values for the charge per DNA 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 DNA 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 r e Verse) that the DNA escapes into the cis chamber and indicates that at low voltages the DNA 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 DNA escape and assume that the DNA 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 DNA does not undergo significant dehybridization prior to escaping into the cis chamber. Under these assumptions the maximal free energy is near state (-(NA - NL),0). 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 p ± reverse k e-(3AG(i,j\-(NA-NL)fl) kre-0AG(i,j\-{NA-NL)fi) + kf6-pAG(iJ\Q,N) = (1 + ^Le(-P&G(-(NA-NL),O\0,M))yl ^ where is an arbitrary ini t ial state, kr is the effective attempt frequency for reverse escape, and kf is the effective attempt frequency for forward escape. Equation 6.16 has a sharp sigmoidal shape at the experimental temperature (288K). Solving Eq . 6.16 for the voltage where P r e Verse = 1/2 yields G + (2N - Nv)svT - \kBT In {^^f) - A * T In ( £ ) V l / 2 ~ \o\(NA-NL) ' ( 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 tran-sitions, we set kr = Li, kf = u>. In the case where the entropic penalty (sv) for storing bases in the vestibule region is <C kB, the numerator of E q . 6.17 is approximately equal to G — kBTln(kf/kr). The difference in the effective attempt frequencies then acts as a correction to the hybridization free energy. If the ratio kf/kr becomes comparable to e130 then Vi/2 is approximately zero, and there is no transition between reverse and forward dominated escape. Equat ion 6.17 suggests that the cross-over between reverse escape and forward escape wi l l occur when the increase in 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), Eq . 6.17 predicts transition voltages of 22mV for the defect free strand, and 14mV for the defective strand. These values are in reasonable agreement with the predictions of the full stochastic model, which predicts 19mV for the defect free strand and l l m V for the defect-bearing strand. When the 2-3 free energy surface is used in Eq . 6.17 the predicted transition voltages for the defect-free and defective strand are 33mV and 29mV, while the full stochastic model gives 23mV and 18mV respectively. Note that Eq . 6.17 used wi th the 2-3 free energy surface gives higher transition voltages than Chapter 6. First Passage Times of Driven DNA Hairpin Unzipping 90 (a) IOOOOF 1000 k-(b) i 0.8 | 0.6 iC 0.4 0.2 0. 0 ( C ) o.( 2 0.4 0.2 0.3 Driving Voltage (V) 0.01 — Defect Free (MFOLD) — - Defect Bearing (MFOLD) •-• Defect Free (2-3) — Defect Bearing (2-3) Drr 0.02 •ing Voltage (V) I • r 0.03 0.04 • V=0.05 • V=0.15 • V=0.30 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 DNA 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 DNA strands escaping into the cis chamber, (c) The probability of having nQ 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. 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 DNA escape. In the exper-iments of Nakane et al, a large Avidin anchor protein was attached to the end of the poly-A tail of the DNA strand, preventing reverse escape. This is equivalent to making the n 0 = — NA states in our model reflecting rather than absorbing. Making this change forces the DNA to dehybridize and escape into the trans chamber, producing an expo-nential 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 DNA segment with a poly-A tail in a a-hemolysin pore. Similarly, Hendrickson et al. [37] found that to hold a completely single stranded DNA 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 DNA hairpin at zero applied voltage occurs on a time scale similar to the time for escape of single stranded DNA 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. C h a p t e r 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 U n z i p p i n g 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 DNA strand and the membrane. This can be seen in Fig. 6.3(c) where -Piast> the probability 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 e r When n 0 < 2N + N\oop so that the transmembrane region is always filled with DNA, sliding of of the DNA is then described by the following master equation: = ~ p ( e - ^ v + l)P(n0,t) + p e - ^ (6.18) where a is the distance moved by the DNA strand in a single sliding event, and A is the charge density along the DNA strand. In the limit of large applied voltage where exp(—/3|A|aV) -C 1, Eq. 6.18 becomes = p ( P ( n Q - l , t ) - P ( n 0 , t ) ) . (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 DNA 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 DNA 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 4 2 , while the 2-3 free energy surface gives a ratio of e 1 8 . 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 4 9 and e2, respectively). We ascribe this difference to the coupling between the dehybridization and the DNA 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 .17 suggests this will occur when the total hybridization free energy is comparable to ksTlnikf/kr). Using kf = u and kr = p and parameters from the fit to the 2-3 free energy, this suggests that the transition between forward and reverse escape would disappear when G « l l /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 DNA 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 NA changes the location of the peak in Chapter 6. First Passage Times of Driven DNA Hairpin Unzipping 94 the mean first passage time curve. Increasing the value of NA means that on average more bases wi l l be located in the high voltage trans chamber in the ini t ia l state of the hairpin. As a result, there wi l l typically be a larger barrier to reverse sliding, the probability of escaping in reverse at any given voltage wi l l be reduced, and the voltage of the peak in the mean escape time wi l l 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 NA up to 60 using the defect-free 2-3 free energy surface and attempt frequencies. Varying has an effect similar to varying NA. When Ni is reduced the number of D N A bases fully inserted into the trans chamber in the ini t ia l state of the D N A typically increases. This results in a larger barrier to reverse sliding, and a reduction in 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 in the vestibule region, results in a larger barrier to reverse sliding as the entropic penalty for reverse sliding increases wi th Nv. However because the fitted values of sv were relatively small, scaling Nv has a minor effect on our results. The qualitative behaviour of the mean first passage time curve remained the same for values of up to 30. However when sv is comparable to ksT, increasing Nv inhibits reverse sliding due to the increased entropic cost of filling the vestibule region. Alter ing Af l o o p , the 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 / 2 . As this quantity 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, in order to reduce the number of fitting parameters in 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 maximum mean first passage times by more than 3 orders of magnitude compared to the results displayed in 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 UJR becomes comparable to p there is an increase of a few millivolts in V i / 2 , and an increase of about half an order of magnitude in the 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 DNA escape previously used to study the DNA translocation process, and set u>R = 0. This could easily be revised in light of new experimental results. 6.4 Summary We have analyzed the processes of translocation of DNA hairpins through narrow trans-membrane 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 DNA 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 DNA sequence on the escape time is minimal. At voltages greater than approximately 25mV the DNA 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 dehy-bridize 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 DNA strands examined ex-hibited 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 DNA 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 transmem-brane 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 bi-ological 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 sim-ulations 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 Model and Simulation Methodology 7.2.1 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 in Ion bearing Nanopores 99 Species cr(nm) e(kJ/mol) 9(e) S P C / E 0.3169 0.6502 q0 = -0.8476, qH = 0.4238 C 0.3283 0.3891 0 N a + 0.2876 0.5216 1.0 K + 0.3250 0.5216 1.0 C a 2 + 0.3019 0.5216 2.0 c r 0.3785 0.5216 -1.0 0.3143 0.6998 -1.0 Table 7.1: The 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]. The 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. The water molecule interacts with the surrounding environment v ia the Coulombic interactions associated with the point charges, and through a Lennard-Jones potential [105] centred on the oxygen atom. Both the partial charges and the Lennard-Jones parameters were originally determined by fit-ting the results of isobaric molecular dynamics simulations to experimental values for the bulk density and dielectric constant of liquid water at 298K [106]. The values for these parameters are listed in Table 7.1. The model used for the ions consists of a point charge that interacts v ia a Coulomb potential and a concentric Lennard-Jones interaction. The parameters for the various ionic species used in this study are also given in Table 7.1. The geometry of the confining cylindrical surface is displayed in F ig . 7.1(a) where the symmetry axis of the cylinder is taken to be the z-axis. The cylinder is modeled as featureless smooth wall using one of two potentials. Following Lynden-Bel l [60] we employ a purely repulsive exponential potential 0waii = 0 o e - ; ( K - r ) , (7.1) to model strictly hydrophobic tubes. In Eq . 7.1, R is the radius of the confining cylinder, and r = \ J X 2 + V2 LS the radial position of the particle. Note that as there is no dependence on z in Eq . 7.1, the confining cylinder is homogeneous in the z direction. The 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] ?wal /OO /»7T / / OO J — 7T \ O 2 2 = eno 7r 63 32' F -9 /2 , -9 /2 ,1 , 3 F ( - 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 Lennard-Jones 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 — 40nm - 2 , the approximate 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 c o n f. 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^altj, VlWh ^IWl) + i=l nw nw / \ 6 / / \ 6 i=l j=i+i \\'o >o\j \\v o ro I sr^ sr^ Qm sr~^ qn i=l j=i+lme{0,Hi,H2} ne{0,HuH2} I'm ' n I nw E E -«=i me{0,HltH2} 47re 0 | r - / - r . m C h a p t e r 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 102 Z = 0 — R = 0.44nm R = 0.54nm - - R = 0.64nm R = 0.74nm T—i—i—nn—i—i—I—i—i—i—i—i—m—1~ / I < I \ I I I _ L J _ o.i 0.2 0.3 r (nm) 0.4 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 kBT is taken to equal the physical radius of the pore. Chapter 7. Water Adsorption in Ion bearing Nanopores 103 In Eqs. 7.3 and 7.4, quantities with the subscripts O, Hi, H2 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, ewwi), (OIW,£IW), and (owwi £ww)y 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 min-imum 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 po-larization 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 TASEP discussed in Chapters 2 through 5, where the steady-state 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 in Ion bearing Nanopores 105 0 0 poo p2n pir p2rr pL p2-n pR ZG = V (nwl)-1^)-^h~3 / . . . / / / . . . / / / e-MH-,*wnw)x „ n -^00 JO Jo Jo J-L Jo Jo n w = 0 n w d3Pld3n Y[ d^dpMdp^dnid3^ i=l (7.5) After integrating out the momentum terms, the partition function can be written n\v=0 Y\ s m ( 8 i ) d Q i d : i r i d 3 r I , nw=Q (2Awe1e2e-i)nw Jo JO JO J-LJO JO n w 1=1 . (7.6) where Ucon{ is the configurational potential energy of the system, j3 = 1/ksT, 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 1 ^/2iimikBT) A = h \j2irniwksT' h ©* = , yJ2ixIikBT 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 ith 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 7. Water Adsorption in Ion bearing Nanopores 106 Peq(nw, {ji}, {f2/},f/) = Z Q 1 A 7 3 ( e - ^ 2 A 3 ^ e i e 2 e 3 ) - n M / e - / 3 t / p ° t ( r i M / ' { r l } ' { ^ } ^ ) x J^[(i3risin(6li)(ir2i(i3r/. i=l (7.7) While Eq. 7.7 in principle provides all the information needed to compute any equilib-rium 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 par-tition function. This is directly analogous to the situation described in Chapter 2 where the large state space of the TASEP made direct solution of the occupancy equations im-practical. 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 Pn+1 = MPn, (7.8) 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, Mtj = 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 Nanopores 107 probability that the transition is successfully completed. 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 wi l l 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. The attempt probability is 5 sm(6)dxdydzdud6d(p °V = wv~ » where 8 is an arbitrary value in the range [0,0.5). The acceptance probability is then Aj = e ^ m a x ( { c ^ 0 ) 1 where 47r 2V \ (7-12) / TT2   Cij = (Ui - Uj) - - kBT\n . \A-wQiQ2®3(nw(]) + i ) y In Eq . 7.12 (Ui — Uj) is the difference in 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/nw, (7.13) and Atj = e - /? m a x (K;>°}) where , m T r , M . rp, f^w^2e3nw(j)\ (7-14) dij = (Ui -Uj)+ pw - kBTIn I J . 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 / 1-25 \ {sm(6)dujd6d(f)\ Aij = e-P™^m-u3)fi}). (7.16) 4. Water Translation: A water molecule currently inside the channel is selected at ran-dom 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 (nw + l ) - 1 - If 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 pro-cess 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 config-uration. 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 M3j = 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, 109 transition attempts were run before 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 109 transition attempts. 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, 0.64 and 0.74nm. 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 kBT 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 Results and Discussion 7.3.1 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 satu-rated water vapour at T = 298A" (pw —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 pw. 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 non-monotonic 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 agree-ment with previous studies [96, 97, 100], the pure-water adsorption-desorption curves for Chapter 7. Water Adsorption in Ion bearing Nanopores 111 -62 -61 -60 -59 -58 -57 -56 -55 -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.w is visibly reduced, particularly in larger 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/cm3), the 0.44nm pore would have approx-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 C h a p t e r 7. Water 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 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/cm3 at 298K), while the ver-tical line gives the chemical potential of saturated water vapour at 298K ( « -58.5kJ/mol). 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 OHO 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. C h a p t e r 7. Water 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 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 dis-played in Fig. 7.5, where the radial distribution function is given by - r x £ E 5 3 (fi ~ f i _ ^) e -^ p o t ( r " / , i ? 1 - -^)d 3 r / fj d 3r f c sin ekdQkd3TZ. i=l j = ljy£i 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 ith 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. An 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 + and C l _ at different values of pw along the filling curve. 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 ion-free, 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 con-taining pure water, and pores containing single ions (pw = —56.0kJ/mol). The ions produce deviations from the pure water structure that are most pro-nounced 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. C h a p t e r 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 118 (a) (b) (c) (d) (e) (f) Figure 7.6: Filling in the 0.44nm pore with a C P ion. (a) p w — —66kJ/mol (nw = 4) (b) p w = —64kJ/mol (nw = 6) (c) the low density state at p w = —62kJ/mol ( n w — 7) (d) the high density state at p w = —62kJ/mol ( n w = 11) (e) p w = —59kJ/mol (nw = 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. 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 ~ 10 _ 6 mol/m 3 , will completely empty the pore of water. As 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 carry-ing a C P and a K + ion, respectively. Note that in all cases there is a strong tendency for C h a p t e r 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 119 (•••••I ! • • • • • * ' •••••I i< ' • • • • • H I !••••••> • • • • • I I • • •P I ' '•••*•» •rr*«i •••••i • • • • • • fe • • • • • I I • • • • • I I '•••••III ( • • • • • I I I • • • • • • I I I '•••••I '•••••I >•••••!! ' • • • • • I I • ••••• • • ? * " ) • ••••••• • • • • • • I • • • • • ! •••••i •••••i • • • • • • !•••••!• !•••••!• 88&1 • * ' J M H j • • £ • • » • M * J « I < •I ' •••••li ' • • • - - • • I : : : : : : : • • • • • • I • • • • • i i • • • • • • I • • • • • i i • • • • • I I >•••••!• _ J f M l i ' • • W « « H '•••> ,)|l< ' • • • • » • • • ( i l P i • • • • • I I was, • ••"»••!« • • ^ • • I l ( • • ^ • • l • • • • • I I • • • • • I I «•••••••< ( • • • • • I I I M I I I I I •n»«i •-•••i • • f c O M •••Mi • ^ * o * •!•••! p ••'•••« !••*••» !•••*••!< > • * • ( " • •••W Hl .<r| • • * H * . ri • • • • • • • • • i i € ^ $?) (s) (§) (J) Figure 7.7: Filling in the 0.44nm pore with a K + ion. (a) p w = -66kJ/mol ( n w = 5) (b) p w = -64.0kJ/mol (nw = 7)(c) p w = -62.0kJ/mol (nw = 9) (d) P w = —61.0kJ/mol ( n w = 10) (e) p w — —59.0kJ/mol (nw = 13) (f) p w = -56kJ/mol (nw = 14). Chapter 7. Water Adsorption in Ion bearing Nanopores 120 (a) (b) (c) (d) (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 w = —64kJ/mol (nw = 8) (c) the low density state at p w = —61.1kJ/mol ( n w = 15) (d) the high density state at p w = —61.1kJ/mol ( n w = 31) (e) P w = —59kJ/mol ( n w — 34) (f) p w = —56kJ/mol ( n w = 40). 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 p o r e s 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 projec-tions 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, C h a p t e r 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 121 Figure 7.9: Filling in the 0.54nm pore with a K + ion. (a) pw = — 66kJ/mol (nw = 7) (b) Pw = —64kJ/mol (nw = 9)(c) the low density state at pw = — 60.2kJ/mol (nw — 17) (d) the high density state at pw = —60.2kJ/mol (nw = 26) (e) Pw = —58.9kJ/mol (nw = 31) (f) pw = —56kJ/mol (nw = 37). an 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 2 + tends to produce higher water densities than N a + , and N a + 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 C h a p t e r 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 122 (d) (e) (f) Figure 7.10: Filling in the 0.64nm pore with a C P ion. (a) p w = -66kJ/mol (nw = 6) (b) P w = —64kJ/mol ( n w — 7) (c) the low density state at p w = —60.6kJ/mol (nw = 19) (d) the high density state at p w = —60.6kJ/mol ( n w = 50) (e) P w = —59.8kJ/mol ( n w = 53) (f) p w = —56kJ/mol (nw = 63). Chapter 7. Water Adsorption in Ion bearing Nanopores 123 (d) (e) (f) Figure 7.11: Filling in the 0.64nm pore with a K + ion. (a) pw = -66kJ/mol (nw = 6) (b) pw = —61.1kJ/mol (nw = 14)(c) the low density state at pw — -59.9kJ/mol (nw = 22) (d) the high density state at pw = -59.9kJ/mol (nw = 49) (e) pw = -58.9kJ/mol (nw = 52) (f) pw = -56kJ/mol (nw = 60). C h a p t e r 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 124 Figure 7.12: Filling in the 0.74nm pore with a CI ion. (a) pw = —66kJ/mol (nw = 5) (b) pw = -64kJ/mol (nw = 6) (c) pw = -60.3kJ/mol (nw = 20) (d) immediately prior to filling pw = —59.9kJ/mol (nw = 25) (e) immediately after filling pw = -59.8kJ/mol (nw = 83) (f) pw — -56kJ/rnol (nw = 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. Chapter 7. Water Adsorption in Ion bearing Nanopores 125 (d) (e) (f) Figure 7.13: Filling in the 0.74nm tube with a K + ion. (a) pw = -66kJ/mol (nw = 6) (b) pw = -64kJ/mol (nw = 7)(c) immediately prior to rilling pw = -59.7kJ/mol (nw = 24) (d) immediately after filling pw = -59.4kJ/mol (nw = 82) (e) pw = -57.5kJ/mol (nw = 88) (f) pw = -56kJ/mol (nw = 92). In (f) a second inner layer of water is beginning to form. C h a p t e r 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 126 10 20 30 40 50 60 15 30 45 60 75 90 Figure 7.14: The total water-ion interaction energy (closed points), and the total particle-wall 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 con-figurational 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 distribu-tions 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 130 r~25 (Ii - 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-water-oxygen peak occurs at approximately 0.27nm,which is nearly identical to the mean water-oxygen-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 nonethe-less 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 (nw) — 36.4), compared to 2.4 bonds per water molecule in Chapter 7. Water Adsorption in Ion bearing Nanopores 131 the K + pores (at pw = — 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 inter-action 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 Nanopores 132 drops with increasing R, and in the limit of R 3> o I W 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 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 = 2ITR2L where R is the geometric radius of the pore. In a large system, far from the boundaries we would expect PL and 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 Cl~ 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 (7.23) E E ( ^ ) + E ^ i ( / 2 - r 0 > i j>i i and the axial pressure given by (7.24) C h a p t e r 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 133 4000 3000 I 2000 • • R = 0.44nm — • R = 0.54nm — R = 0.64nm *—* R = 0.74nm -1—I—I—I—I—I—I—I-Figure 7.19: (a) The radial pressures in K + (closed symbols) and Cl~ (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 ~ 102atm. In the case of PL, the ideal term is typically larger than the observed pressure, 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 w corresponding to the filling regions of the isotherms in Fig. 7.3. These 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 approx-imately 8 x 10~3 atm and compare this with the radial pressures of Fig. 7.19(a), it is 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 2GPa [115], and recent 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 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 PR and PL are expected to depend on the radius of the confining 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 PL and PR would be independent of L. With these assumptions and under conditions of fixed p,w, T, and R dfi = -SdT - A-nRLPRdR - 2nR2PLdL - nwdp (7.25) dQ = -2ixR2PLdL. (7.26) As Pi is independent of L we can trivially integrate to find that tt(pw, T, R, L) = -2nR2LPL = -PLV. (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 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 nU(S, V,L,nw) = U(nS,r]V,riL,r]nw) (7.28) PL(L = e)-P P < 1. (7.29) Let QQ = —PV and then consider the quantity (7.30) Now let 7 = max (\PL(x)\) 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). < 1 + P (7.31) 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 expec-tation, the trivial integration of the grand potential is only valid in the thermodynamic limit of large L. Now consider Figure 7.20 which shows PL as a function L in a pure water pore and in a K + 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 PL averaged over these oscillations is close to the value of PL at large L. Conversely, the ion bearing pore shows a large nearly monotonic decrease in PL as L is increased. The results of Fig. 7.20 suggest that for a K + 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 po-tentials 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 In Eq. 7.32, pw 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(pwl\v, L, T) = — kBTm(Zo). In the discussion that follows, pff = — 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 po-tentials, and demonstrates that in pores with R > 0.54nm the per particle entropy is (7.32) C h a p t e r 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 N a n o p o r e s 138 -100 -125 -150 -175 -200 E -63 -62 -61 -60 -59 -58 " " " -63 -62 -61 -60 -59 -58 G -150 -200 -250 -300 ~ 3 5 - 6 3 -62 -61 -60 -59 -58 -63 -62 -61 -60 -59 -58 H w (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 | 6 8 10 12 14 10 15 20 25 30 35 I 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 ~ R2. 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 water-water interaction can be justified by noting that the water molecules can be approximated as dipoles with moments of approximately 7.8 x 10~ 3 0C • m. Setting r to be the distance between two water molecules, the water-water interaction energy would then scale as r~ 3 if the relative orientation of the molecules were fixed, and as r~ 6 if there is thermal averaging over the orientations. If the interfacial energy goes as ~ R2, then if R is fixed and L is scaled, the interfacial energy will decrease as a fraction of the total energy in the system at least as fast as ~ L~l, assuming the total energy in the filled state is proportional to the total volume. Based on these scaling arguments, we expect that as L is increased holding R fixed C h a p t e r 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 140 80 60 40 20 0 •— L= 1.50 nm * — L = 3.10nm :-. ( a ) 1 i . . . > . . , i . . . . i . . , . ) , . . , r --66 -64 -62 -60 -58 -56 -54 0 0.1 0.2 0 0.1 0.2 r (nm) o i 1 1 1 1 1 (d) Hydrogen A • 1111111 0 0.1 0.2 0.3 0.4 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 (pw = —58.2 kJ/mol). 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 water-water 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 C h a p t e r 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 141 K 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 c a p = 2^0e-'Lcosh(/^). (7.33) 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. An 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 ad-sorption 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 chlo-ride 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 103atm. 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 examina-tion 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 TASEP showed that the TASEP with extended particles retained the three-phase character of the standard TASEP. By explicitly enu-merating 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 mo-tion [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 TASEP lat-tice. 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({ni},L) ( E i n i ) ! / i V - E i " i ( r f i - l ) ) (8.1) Chapter 8. Conclusions and Suggestions for Future Work 147 where 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 where Oj is the bulk density of particles of size j. The current J could then be maximized with respect to the ox given specific constraints such as Oi/oi+1 = ai/ai+i in the case where am < j3n for all m and n. This preliminary analysis suggests that there are no 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 p2- 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 TASEP 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 TASEP 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 p2 was reduced, and a loss of symmetry 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 TASEP could be performed, both for intrinsic interest in the behaviour of this variant of the TASEP 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. An interesting extension would be to consider the case where the spatial variation in the hopping rates was smooth, with a gradual change in (8.2) 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 TASEP 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 TASEP 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 re-duction 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 signifi-cant 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 DNA through nanopores is another biologi-cally related transport process, and was investigated in Chapter 6. Following recent exper-imental work, we considered the voltage driven translocation of DNA 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 DNA hairpin was found to escape the a-hemolysin channel by sliding in a direction Chapter 8. Conclusions and Suggestions for Future Work 149 opposite to the one favoured by the applied voltage. 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 compa-rable 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 DNA seg-ment 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 DNA 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 DNA 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 103 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 anion-cation 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' G u i d e : Solu-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 w i t h I m p l i c i t y l y R e s t a r t e d A r n o l d i Methods (SIAM, Philadelphia, Pennsylvania, 1998). [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 185, 511 (1996). [9] T. J. Minehardt, N. Marzari, R. Cooke, E . Pate, P. A. Kollman, and R. Car, Biophysical Journal 82, 660 (2002). [10] J. D. Lawson, E . Pate, I. Rayment, and R. G. Yount, Biophysical Journal 86, 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 69, 667 (1992) . [14] S. A. Janowsky and J. L. Lebowitz, Physical Review A 45 , 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 72, 277 (1993). [17] A. B. Kolomeisky, Journal of Physics A 31 , 1153 (1998). Bibliography 153 [18] T. Nagatani, Journal of Physics A 26, 6625 (1993). [19] S. Sandow, Physical Review E 50, 2660 (1994). [20] M. Schreckenberg, A. Schadschneider, K. Nagle, and N. Ito, Physical Review E 51, 2939 (1995). [21] F. H. Essler and V. Rittenberg, Journal of Physics A 29, 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 UK, 1997). [23] G. Tripathy and M . Barma, Physical Review Letters 78, 3039 (1997). [24] T. Sasamoto and M . Wadati, Journal of Physics A 31 , 6057 (1998). [25] G. Tripathy and M . Barma, Physical Review E 58, 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 31 , 6911 (1998). [28] M. Bengrine, A. Benyoussef, H. Ez-Zahraouy, J. Krug, M. loulidi, and F. Mhirech, Journal of Physics A 32, 2527 (1999). [29] V. Karimipour, Physical Review E 59, 205 (1999). [30] V. Popkov and G. M. Schiitz, Europhysics Letters 48 , 257 (1999). [31] K. M. Kolwankar and A. Punnoose, Physical Review E 61 , 2453 (2000). [32] Z. Nagy, C. Appert, and L. Santen, Journal of Statistical Physics 109, 623 (2002). [33] T. Chou, Biophysical Journal 85, 755 (2003). [34] L. B. Shaw, R. K. P. Zia, and K. H. Lee, Physical Review E 68, 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 85, 3057 (2000). [38] A. Meller, L. Nivon, and D. Branton, Physical Review Letters 86, 3435 (2001). [39] S. Howorka, S. Cheley, and H. Bayley, Nature Biotechnology 19, 636 (2001). [40] W. Vercoutere, S. Winters-Hill, H. Olsen, D. Deamer, D. Haussler, and M. Akeson, Nature Biotechnology 19, 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. Li , 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 i n S t a t i s t i c a l Physics (Oxford University Press, Oxford, United Kingdom, 1999). [69] L. B. Shaw, J. P. Sethna, and K. H. Lee, Physical Review E 70, 021901 (2004). [70] D. Chowdhury and J. S. Wang, Physical Review E 65, 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 44 , 6375 (1991). [73] M. Bulmer, Genetics 129, 897 (1991). [74] G. Marais and L. Duret, Journal of Molecular Evolution 52, 275 (2001). [75] J. R. Powell and E. N. Moriyama, Proceedings of the National Academy of Sciences of the United States of America 94, 7784 (1997). [76] C. G. Kurland, FEBS 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 (Springer-Verlag, New York, New York, 2003). [78] M. Robinson, R. Lilley, S. Little, J. S. Emtage, and G. Yamamoto, Nucleic Acids Research 12, 6663 (1984). [79] M. A. Sorensen, C. G. Kurland, and S. Pedersen, Journal of Molecular Biology 207, 365 (1989). [80] D. A. Phoenix and E. Korotkov, FEMS Microbiology Letters 155, 63 (1997). [81] W. Konigsberg and G. N. Godson, Proceedings of the National Academy of Sciences of the United States of America 50, 687 (1983). [82] G. F. T. Chen and M. Inouye, Nucleic Acids 18, 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 67, 509 (1998). [85] A. F. Sauer-Budge, J. A. Nyamwanda, D. K. Lubensky, and D. Branton, Physical Review Letters 90, 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 274, 1859 (1996). [88] P. G. de Gennes, Scaling Concepts i n P o l y m e r Physics (Cornell University Press, Ithaca, New York, 1979). [89] A. Hanke and R. Metzler, Journal of Physics A 36, L473 (2003). [90] G. Altan-Bonnet, A. Libchaber, and O. Krichevsky, Physical Review Letters 90, 138101 (2003). [91] J. Chuang, Y. Kantor, and M. Kardar, Physical Review E 65, 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, 18204 (2004). [95] D. Nicholson and N. Quirke, Molecular Simulation 29, 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 21 , 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 100, 1420 (1996). Bibliography 157 [105] M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Oxford University Press, Oxford, United Kingdom, 1987). [106] H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, Journal of Physical Chem-istry 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 Mechanics (Addision-Wesley, Reading, Massachusetts, 1980), 2nd ed. [110] D. A. McQuarrie, Statistical Mechanics (Harper Collins, New York, New York, 1976). [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 Physics (John-Wiley & Sons, New York, New York, 1998), 2nd ed. [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 Monte Carlo Algor i thm 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#fcrfcJ-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 TASEP 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 TASEP 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 A p p e n d i x 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 (Z4) is the selected move. 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 B K L simulation is the so-called binary heap. An example can be viewed in Fig. A . l . In Fig. A . l each circle represents a single node 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 nx, then any nodes on a higher level in the heap that are connected to node are said to be children of m . In a binary heap, a node may have a maximum of two children which we term the left child and the right child. Similarly, a node nx may be connected to a maximum of one node that is at a lower level in the heap, and is called the p a r e n t of m . Nodes that have no children are called leaves of the heap, while the single node that has no parent is called the root or head 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 TASEP 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, / e n try = 0, otherwise /entry = 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 TASEP state. Also note that the number of levels in the binary heap is [\og2(2N + 1)J + 1. 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 TASEP state. First note that for a TASEP with N sites, the associated binary heap contains 2N + 1 floating-point values. For values of N ~ 104, it is convenient to store the entire binary 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 TASEP to escape its current state. Essentially the algorithm then executes a binary search for the transition with a rate r{ such that Efc=orfc < s < J2k=ork- The algorithm then returns the index 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. An 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(log2(Ar)) time, 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 104 sites with extremely good performance. While not absolutely 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 162 Algorithm 1 L C H I L D ( Z ) Requires: An 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 R C H I L D ( Z ) Requires: An 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 P A R E N T ( Z ) Requires: An 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\ Appendix A. An implementation of the Bortz-Kalos-Lebowitz Monte Carlo Algorithm 163 Algorithm 4 L E A F I N D E X ( Z , N) Requires: An integer 1 < i < N corresponding to the index of a lattice site (i) Requires: An integer equal to the length of the TASEP lattice (A'') Returns: The array index of the leaf corresponding to lattice site i 1: return N + i + 1 Algorithm 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 TASEP with the correct movement rate stored for lattice site i 1: I N T E G E R j <— L E A F l N D E X ( i ) 2: Heap[j] <— rate 3: j <— P A R E N T ( J ' ) 4: while j > 1 do 5: Heap[j] <- Heap[LCmLD(j)] + F e a p ( R C H l L D ( j ) ] 6: j <— P A R E N T ( j ) 7: end while 8: return Algorithm 6 S E L E C T T R A N S I T I O N ( A ^ , Heap) Requires: The number of sites in the TASEP 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: while i < k do 5: if s < Heap[LEFTCHILD(i)] then 6: i < — L E F T C H I L D ( Z ) 7: else 8: s <— (s - i / eap [LEFTCHILD ( i ) ] ) 9: i ^ R l G H T C H I L D ( z ) 10: end if 11: end while 12: i <— (i — k) {Injection events return 0} 13: return i Appendix A. An implementation of the Bortz-Kalos-Lebowitz Monte Carlo Algorithm 164 Algorithm 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: An integer equal to the number of sites in the TASEP lattice, (N) Requires: An integer array of length N with each element equal to 0, (Lattice) Requires: A floating-point array of length 2Af + 1 with each element equal to 0, (Heap) Requires: An integer equal to the maximum number of simulation iterations, (M) I N T E G E R i I N T E G E R j F L O A T 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 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

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>
                        
                    
IIIF logo 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

Comment

Related Items