Science, Faculty of
Physics and Astronomy, Department of
DSpace
UBCV
Wong, Wesley Philip
Master of Science - MSc
University of British Columbia
A Monte Carlo simulation was used to investigate the thickness compressibility and shape properties of a tethered polymer loop. Specifically, the dependence on chain length, bending stiffness (i.e. persistence length) and excluded volume interactions was examined. Power law relationships were determined when appropriate. For a flexible phantom chain of 90 segments, reasonably good agreement was found with the predictions of an analytical random flight model based on the diffusion equation. By examining the transition from flexible to semi-flexible, rough bounds were placed on the domain in which a worm-like chain can be well described by an equivalent freely jointed chain. As the bending stiffness was increased but the contour length held fixed, the polymer loop crossed over into the "stiff elastic rod" regime. This was marked by a change in the relationship between force and gap distance, as well as an alteration of shape properties. Finally, the single polymer loop model was compared to a simulation of the red blood cell spectrin network by Boal [1]. Many of the features of network geometry and elastic thickness compliance could be well approximated by the single chain model.
5473803 bytes
application/pdf
Modelling a Tethered Polymer Loop: Shape and thickness compressibility by Wesley Philip Wong B.Sc, The University of British Columbia, 1999 A THESIS SUBMITTED IN PARTIAL F U L F I L M E N T OF T H E REQUIREMENTS FOR T H E D E G R E E OF M A S T E R OF SCIENCE in T H E F A C U L T Y OF G R A D U A T E STUDIES (Department of Physics and Astronomy) We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA September 26, 2001 © Wesley'Philip Wong, 2001 In presenting this thesis in partial fulfilment of the requirements for an ad-vanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Physics and Astronomy The University of British Columbia Vancouver, Canada Abstract ii ABSTRACT A Monte Carlo simulation was used to investigate the thickness compress-ibility and shape properties of a tethered polymer loop. Specifically, the dependence on chain length, bending stiffness (i.e. persistence length) and excluded volume interactions was examined. Power law relationships were determined when appropriate. For a flexible phantom chain of 90 segments, reasonably good agreement was found with the predictions of an analytical random flight model based on the diffusion equation. By examining the transition from flexible to semi-flexible, rough bounds were placed on the domain in which a worm-like chain can be well described by an equivalent freely jointed chain. As the bending stiffness was increased but the contour length held fixed, the polymer loop crossed over into the "stiff elastic rod" regime. This was marked by a change in the relationship between force and gap distance, as well as an alteration of shape properties. Finally, the sin-gle polymer loop model was compared to a simulation of the red blood cell spectrin network by Boal [1]. Many of the features of network geometry and elastic thickness compliance could be well approximated by the single chain model. Contents iii C O N T E N T S Abstract ii List of Tables vi List of Figures vii Acknowledgements 1. Introduction 1 1.1 Comments on Organization of the Thesis . 4 2. Computational Methods 6 2.1 Introduction to the Monte Carlo Method 6 2.1.1 Integration by sampling 7 2.1.2 The Markov Chain 8 2.2 Implementations of the Monte Carlo method and Applica-tions to Physics 11 2.2.1 Statistical physics and ensemble averages 12 2.2.2 The many-body problem 13 2.2.3 Standard Monte Carlo algorithm 14 2.2.4 Smart Monte Carlo algorithm 14 2.3 Generation and Analysis of the Data 16 2.3.1 The Simulation 16 2.3.2 Additional Analysis 17 3. Isolated Chain in Solution 19 3.1 Single Chain Polymer Models 19 3.1.1 The freely jointed chain 19 3.1.2 The worm-like chain 20 3.1.3 Models with excluded volume and self-avoidance . . . 2 1 3.1.4 The computer model 22 3.2 Shapes and dimensions in large N limit 25 Contents iv 3.2.1 Global length scales 26 3.2.2 Shape properties 38 3.3 Persistence length 43 3.3.1 Results and Analysis of Computer Simulations . . . . 44 4. Single Chain Tethered to an Interface 50 4.1 Single Polymer Loop Models 50 4.1.1 The Diffusion Model 51 4.1.2 The computer model 54 4.2 Shape Properties as a Function of N 55 4.2.1 Results and Analysis of Computer Simulations . . . . 55 4.3 Shape Properties as a Function of Bending Stiffness 67 4.3.1 Results and Analysis of Computer Simulations . . . . 68 5. Single Loop under Parallel Confinement 77 5.1 Models and Theory 77 5.1.1 The Diffusion Model 77 5.1.2 Untethered polymer models 80 5.1.3 The computer model 81 5.2 Results and Analysis of Computer Simulations 82 5.2.1 Shape distortions and compressibility of a long flexible loop 82 5.2.2 Dependence of compressibility on N and kc under strong confinement 87 6. Concluding Discussion 93 6.1 Comparison with a Full Network Simulation 93 6.2 Discussion on Elastic Thickness Compressibility 97 Bibliography 103 A. Glossary of Symbols 108 B. Basic Concepts of the Monte Carlo Method 110 B . l Simple sampling 110 B.2 Formalism: Probability spaces, random variables, and distri-bution functions I l l B.2.1 Probability spaces I l l B.2.2 Conditional probability 112 B.2.3 Random variables and distribution functions 112 B.2.4 Expectation values and functional composition . . . . 114 Contents v B.2.5 The law of large numbers 114 B.3 Importance Sampling 115 B.4 The Markov Chain 116 B.4.1 Heuristic argument 117 B.4.2 Algorithm 118 B.4.3 Continuous random variables 119 List of Tables vi LIST O F T A B L E S 3.1 Contour length fit for various models 28 3.2 Scaling laws for (r^e) and (rj?) for various models 30 3.3 Scaling laws for (r1e) and (r^) for chains with kc ^ 0 33 3.4 Mean large N shape properties for various models 39 4.1 Power law fits for single polymer loops 57 4.2 Characterization of the height distribution of single polymer loops 59 6.1 Scaling laws for L/(z) and AC/(A) for tethered polymer net-works 95 6.2 Scaling laws for the dimensionless elastic constant /3YL3 for tethered polymer networks 98 List of Figures vii LIST O F F I G U R E S 1.1 Red blood cell membrane architecture 2 3.1 Mean contour length vs. number of segments 29 3.2 (r2e) and (r2) vs. number of segments 31 3-3 (r1e) prefactors and exponents vs. sphere radius 32 3.4 (r-ge) and (r2) as a function of N for P P C simulations with non-zero kc 34 3.5 (r2e) and (r2g) for PPC100 chain with kc = 6 35 3.6 Prefactor of linear fit to {r2e) vs. N for various settings of kc 36 3.7 (r2ee) and (r2g) as a function of N for SA100 chain with kc = 1 37 3.8 Mean asphericity vs. number of segments 40 3.9 Mean S3 vs. number of segments . 41 3.10 Mean ratios of principle moments of inertia vs. number of segments 42 3.11 Tangent-tangent correlation function 46 3.12 Persistence length as a function of kc for PPC100 47 3.13 Persistence length as a function of kc for SA100 48 4.1 (r-gg) and (r2) vs. number of segments for a polymer loop . . 56 4.2 (ree) and (top height) vs. number of segments for a polymer loop 57 4.3 Variance of top height vs. number of segments for a polymer loop 59 4.4 Top height modulus vs. number of segments for a polymer loop 60 4.5 Height distribution for PPC100 polymer loops 61 4.6 Height distribution for SA100 polymer loops 62 4.7 Descriptive statistics for the height distribution of polymer loops, as a function of segment number 63 4.8 Mean height vs. number of segments for a PPC100 polymer loop 64 4.9 Effective Nb2 vs. number of segments for a polymer loop . . 66 4.10 (r-gg) and (r2) as a function of ked for a polymer loop 69 List of Figures viii 4.11 (ree) and (top height) as a function of keff for a polymer loop 70 4.12 Variance of top height as a function of kefi for a polymer loop 71 4.13 Top height modulus as a function of fceff for a polymer loop . 72 4.14 Height distributions for polymer loops with different values of kc 73 4.15 Descriptive statistics for the height distribution of polymer loops as a function of £;eff 74 4.16 Descriptive statistics for the height distribution of polymer loops as a function of keff 76 5.1 Free energy and force vs. gap distance for the Diffusion Model 79 5.2 Height measurements as a function of force for a confined polymer loop 84 5.3 Height distributions for confined polymer loops 85 5.4 (r^e) and (r^) vs. mean gap distance for a polymer loop . . . 86 5.5 Mean top height and mean gap distance vs. compressive force for confined polymer loops of various lengths 88 5.6 Power law prefactors and exponents of compressive force vs. number of segments for a confined polymer loop 89 5.7 Mean top height and mean gap distance vs. compressive force for confined polymer loops of various bending stiffnesses . . . 90 5.8 Power law prefactors and exponents of compressive force vs. kes for a confined polymer loop 91 5.9 Power law prefactors and exponents of compressive force vs. ked for a confined polymer loop II 92 6.1 Scaled reciprocal thickness and surface density vs. number of segments for tethered polymers in networks 96 6.2 Dimensionless elastic constant / 3 Y L 3 vs. number of segments for tethered polymer networks 98 6.3 Scaled stiffness per chain (3YLA/3 vs. number of segments for tethered polymer networks 99 6.4 Height measurements as a function of force for a confined polymer loop (revisited) .' 101 6.5 Height measurements as a function of stress for a confined polymer loop 102 A cknowledgem en ts ix A C K N O W L E D G E M E N T S I would like to express my deepest appreciation to Dr. Evan Evans, who has been a wonderful teacher and mentor, far above the call of duty. There are many people whose assistance and support have made this work possible, and I am grateful to all of them. I would especially like to thank: Myer Bloom, for getting me started on this research path, and for taking the time to read through this thesis; Volkmar Heinrich, for his patience and help with final revisions; Andrew Leung and Wieslawa Rawicz, for helping to make my time at U B C so enjoyable; My buddies Michael Forbes, Mark McCann, Ari Turner, and Joel Erickson, for many fruitful (or at least entertaining) discussions; My family, Dr. Paul, Dr. Lilian, and Dr. P. Austin Wong; I don't know where I would be without their support; Jesus Christ . . . for everything good. This work was supported by the NSERC of Canada. 1. Introduction 1 1. I N T R O D U C T I O N Polymers play an important role in many areas of science. They are used in nature as the building blocks of life, essential for structure, function and information storage within living cells (e.g. microtubules, proteins, DNA). Within the past century, the synthesis of polymers in the laboratory has re-sulted in many technological innovations (e.g. first synthetic fibre, Chardon-net, 1887; first plastics, Baekeland, 1909) [2]. The study of polymers is an active field in both basic and applied research. One system that has garnered much attention is the membrane cytoskele-ton of the red blood cell (RBC). The cytoskeleton is a network of protein filaments that is attached to the cytoplasmic (interior) side of the plasma membrane (the encapsulating interface that acts as the cell's "skin"). It plays a vital role in the mechanical and functional properties of the cell. The human R B C membrane is a good model system for studying a cytoskele-ton, because its design is prototypical of the encapsulating membranes of eukaryotic cells. The cytoskeleton in the R B C consists primarily of a tri-angulated network of spectrin molecules. These molecules form tetramers with a contour length of approximately 200 nm, and a persistence length1 of approximately 10 nm [3-5]. A depiction of the red cell membrane ar-chitecure is presented in Figure 1.1. While a great deal is known about the R B C membrane (as per the review article [3] or the molecular biology text-book [6]), questions remain as to how the pieces interact, particularly in the dimension perpendicular to the lipid bilayer. Recently, experiments of the thickness compressibility of the membrane have been conducted [7] using an ultra-sensitive force probe [8]. In order to understand the contribution of the spectrin network to the thickness and elasticity of the R B C membrane, additional research is required on the physics of tethered polymer networks. The development of models for describing polymers began in 1920 with Staudinger, who introduced the idea that polymers are long chains of re-peating units [2]. Foundational work in the physics of polymers was es-1 Persistence length is a measure of the flexibility of a filament, quant i fy ing the dis-tance a long the molecule between uncorre lated tangents. T h i s wil l be discussed further in Chap t e r 3. 1. Introduction 2 F F i g . 1.1: Schematic representation of red cell membrane architecture reproduced from Picar t and Discher (1999) [7, 9] that depicts the l ip id bilayer (thin-lined slab) studded wi th a superficial forest of glycoproteins and supported by a subsurface network of spectrin tetramers which is linked to the bilayer by actin/4.1 junctional complexes. The complex of band 3 /ankyr in (light gray) is sketched as putatively bound to spectrin between the junctional complexes. Also shown is hypothetical confinement of a spectrin loop by a rigid wall . The wall is pushed by a spring designed to report the steric repulsion from the loop. The geometric and elastic properties of the tethered loop have been the major focus of the work reported in this thesis. 1. Introduction 3 tablished by Flory, Debye, Kuhn and Kramers in the 30s and 40s [10-12]. S. F. Edwards pioneered the application of the modern theoretical tech-niques of many-body theory to polymer physics in the 70s [13]. More re-cently, theoretical interest was revitalized by de Gennes and des Cloiseaux, who established the relationship between polymer physics and critical phe-nomena [10, 14, 15] (e.g. the polymer-magnet analogy). An observation that has facilitated much theoretical work in this area is that in many cases the macroscopic properties of long chain molecules depend only on a few key parameters, but are otherwise independent of the microscopic details. Thus, universal properties of polymers can be discovered by studying the behaviour of very simple models. For example, in the extreme case of long flexible polymers, diffusion/random walk theory is an effective modelling tool [7, 16-18]. However, these models do not apply to semi-flexible poly-mers, in which the persistence length and contour length are comparable. Far fewer analytical calculations have been established in this regime, where both entropic and energetic factors are significant. Moreover, even less work has been done on polymer chains in networks tethered to a surface. Analytical models have been constructed in order to understand the vis-coelastic properties of semi-flexible polymer networks. Two such models are the "affine" [19, 20] and "tube" [21] models. The affine models are based on the assumption that under material extension, the polymer filaments un-dergo an affine deformation on length scales greater than the entanglement length.2 In other words, the deformation of the polymers on a mesoscopic scale mirrors the deformation on the macroscopic scale. This assumption al-lows macroscopic properties to be predicted from the longitudinal response of a polymer filament. On the other hand, the "tube" model is based on the effective confinement of polymer filaments that arises from network entan-glements and other chain-chain interactions [21]. Both models use a single filament picture in order to describe bulk properties. Two issues that are difficult to address analytically are the excluded vol-ume effect, and the tethering of polymers to an interface. The excluded vol-ume effect refers to the swelling of a polymer due to repulsive self-interactions (e.g. sterics), which typifies real polymers in good solvent. The effect of tethering a polymer to an interface is important for the R B C cytoskeleton network, as well as other systems in nature. An understanding of such sys-tems also has industrial applications, for example in steric stabilization of colloids in paint [22, 23]. 2 T h e entanglement length describes the average distance between steric constraints or cross-links. 1. Introduction 4 A detailed method for studying steric interactions and tethered archi-tectures is through computer simulations. For example, these techniques were recently used to simulate the R B C spectrin network [1, 24-26]. The standard approach is to model the polymer (e.g. a spectrin tetramer) as a string of beads connected by string-like tethers (e.g. [1]), or to course grain the molecule even further (e.g. [26]). While random flight-analytical models are effective for studying long flexible chains, it is not clear how well these models work for the spectrin network of the R B C . To my knowledge, there have only been a couple of computer simulations that have examined the elasticity and geometry of a tethered polymer network [1, 25]. Estimates of a transverse modulus were obtained from these simulations, based upon the mean and variance of the bead heights, but the network was not examined under conditions of parallel confinement. Attempts to compare the simu-lation results with experiments on the R B C thickness compliance appeared to yield a modulus much greater than experimentally measured [7]. Thus, as depicted in Figure 1.1, the main objective of this thesis has been to to characterize the shape and mechanical properties of a polymer tethered to an interface. 1.1 Comments on Organization of the Thesis In this work, I have investigated how the shape dimensions and elastic thick-ness compressibility depend on the parameters of bending stiffness, excluded volume and length. Power law exponents are given when appropriate. In addition, some characterization of the domain of applicability of computa-tional random flight models and diffusion models is provided. The purpose of this study has been to provide a. first step towards developing a predic-tive model of the shape and thermodynamic compliance of tethered polymer networks. In the dilute regime, where interactions between chains are weak, the single chain results may be used to approximate a network. Further-more, when building computational network models, it is helpful to have a well-characterized building block to work from. This can aid in deciding whether an explicit bending stiffness potential is required. The next chapter describes the computational methods that were used in this investigation. A presentation of the Monte Carlo method is included. The following three chapters reflect the fact that the model was built in stages. The three stages are: the isolated chain in solution, the tethered polymer loop, and confinement of the tethered polymer loop. Analytical predictions are presented and compared to computational results in each of these chapters. The thesis ends with a concluding discussion. Appendices 1. Introduction 5 include an expanded introduction to the Monte Carlo method, and a table of frequently used symbols and abbreviations. 2. Computational Methods 6 2. C O M P U T A T I O N A L M E T H O D S In this study, the properties of polymer systems were investigated through the use of computer simulations. In order to assess the reasonableness and validity of the simulations, the results were compared to established analyt-ical models whenever possible. The strategy was to build up the model in stages, from simple systems which could be checked against exact results, to more complex systems of confined, tethered polymers. There are many different methods for simulating physical systems. As I was primarily interested in investigating static equilibrium properties, I used a Monte Carlo method instead of standard Brownian dynamics techniques. However, the Smart Monte Carlo technique which was used for many of the simulations, can also be used to calculate the time evolution of ensem-ble averages with an explicit time variable. It is in some sense a hybrid technique. In the first section of this chapter, the basic ideas and theory behind the Monte Carlo method are presented. An expanded version of this sec-tion, which includes more mathematical details, is provided in Appendix A. The second section presents various implementations of the Monte Carlo technique, including the Smart Monte Carlo algorithm. Also included is a physical interpretation of the prescription. Finally, the simulation which was used in the present study is discussed, as are the tools used to generate and analyze the simulation data. 2.1 Introduction to the Monte Carlo Method The Monte Carlo Method refers to a collection of stochastic techniques de-veloped by Metropolis, Ulam, and von Neumann in the 1940s. For this reason it is sometimes called the Metropolis method. One of the central problems that this method addressed was the numerical calculation of inte-grals through the use of statistical sampling. This is the technique that will be presented here. The discussion is based primarily on references [27], [28] and [29]. It will be assumed that the reader is familiar with standard proba-bilistic concepts such as random variables (r.v.s), and probability densities. 2. Computational Methods 7 For a review of the relevant aspects of probability theory, see Appendix A section B.2. 2.1.1 Integration by sampling Consider the problem of approximating the 1-dimensional integral: rd f(x)dx (2.1) o If this integral is well-defined, it can be approximated in the standard way by selecting N equally spaced points in the interval [0,d], then evaluating the function / at these locations. This yields the expression, , s i l i / f r o aPP \T ^ (2-2) where S = {xt} is the set of uniform x coordinates. Clearly, as N approaches O O ) lapp approaches / . The basic idea behind the Metropolis approach is the following: the X j ' s do not have to be chosen in a grid in order to yield a meaningful approxima-tion; rather, they can be chosen randomly between 0 and d. Then if the new set [xi] is selected with uniform probability within this interval, equation (2.2) still represents a valid approximation of the integral I as N becomes large. In the language of probability, {xi} is said to contain N instances of the random variable X, where X has probability density p(x) = 1/d on the interval [0,d] and zero probability density elsewhere. Furthermore, Iapp/d is an approximation of the mean value {f(X)) = I/d. This approximation method, in which p(x) is constant, is known as simple sampling. Now suppose that p(x) is any integrable function. Then the mean value of f(X) is given by: /+oo f(x)p(x)dx (2.3) •oo So in order to approximate / , let us require p(x) to be non-zero on the domain of integration [0, d], and zero elsewhere. In addition, let us define the function: 2. Computational Methods 8 Then the average value of the random variable g(X) yields the desired in-tegral: As before, this integral can be approximated by averaging specific in-stances of g{X), with the probability distribution of X given by p(x). Notice that dividing by p(x) takes care of the normalization. When the probability density of X is not constant over the domain of integration, the technique is called importance sampling. This method can be more efficient than sim-ple sampling as p(x) can be chosen to preferentially sample regions of the domain of integration where f(x) is the greatest. It is clear that this technique is not limited to integration over one di-mensional intervals, but can be applied to much more general domains. Here is a summary of the basic strategy: to integrate a function / over some do-main, choose a random variable X with probability density p non-zero on this domain, and zero elsewhere. Then the integral is the average value of (f/p)(X), and this average can be approximated by sampling the random variable. In order to carry out the Monte Carlo method, it is necessary to construct a random variable with a specified probability density. A powerful technique for accomplishing this was developed by Metropolis et al. The basic idea is to form a sequence of random variables, specifically a Markov Chain, which converges to the desired r.v. In this subsection, I will present a brief introduction to Markov chains, followed by a description of the prescription. Consider a sequence of random variables indexed by a "time variable": {Xi : i G Z>o}. To simplify the present discussion, it will be assumed that the index is discrete. Furthermore, let us also assume that the random variables are discrete, and map to the (ordered) state space S. Let the probability mass function of each variable Xi be denoted by pi. If each r.v. depends on the previous variables in a well-defined manner this sequence is called a random process. If in addition, the (conditional) distribution of +00 (2.5) 2.1.2 The Markov Chain 2. Computational Methods 9 Xi depends only on the random process is called a Markov Chain. Formally, this requirement is 1 P(Xn = s\X0 = x0,-.., Xn-i = xn-.i) = P(Xn = s\Xn-i = xn-\) Vn > 1 and all s, XQ, . . . , xn-\ £ S (2-6) Now suppose we want to construct a random variable X with probability mass function p. The prescription given by Metropolis et al. yields a Markov Chain such that for sufficiently large N, Xjy has a distribution which is sufficiently close to p. The key to this construction is the condition of detailed balance: p(i)Wi^f = p(f)wf^, V i j e s (2.7) where Wa->b represents the transition probability from state a to state b, that is Wa^b = P(Xk = b\Xk^ =a)Vk (2.8) When this condition is satisfied, we would like to have the following conver-gence: for all e > 0, there exists an N > 0 such that ||P(-X"n = a;) — P(X = x)\\ < e, whenever n > N. Although there are some subtle issues which will not be addressed here, for practical purposes this convergence can usually be taken to hold. A simple heuristic argument to show that this is reasonable is presented in Appendix A, section B.4.1. Algorithm In practice, one can generate an instance of the desired Markov Chain in the following way: 1. Initialization. Choose a starting state XQ = a, where Xi refers to the instance of the r.v. Xi. Also, initialize the index: i = 0. 2. Make a trial move. Consider the trial move as an instance of the ran-dom variable -Xt r ; a i . The conditional probability of choosing X t r i a i = / given Xi = X j is denoted by P{Xu-m[ = / | X c u r r e n t = x^. We write - ^ c u r r e n t instead of Xi to emphasize the fact that this probability is usually chosen so that it does not depend explicitly on i. 1 The conditional probability that event A occurs given that event B occurs is denoted by P(A\B) = P(A n B)/P(B). 2. Computational Methods 10 3. Choose or reject this move according to some prescribed probability, P(Accept move? | A ^ c u r r e n t = & Atrial = /)• 2 As before, we are assuming that this probability does not depend explicitly on the index i. If the move is accepted, let Xi+\ = / . Otherwise, let Xi+i = x\. 4. Increase i by one and go back to step 2. This procedure gives us a lot of freedom both in choosing a trial move, and in setting the probability of acceptance. The primary condition that must be satisfied is detailed balance (equation 2.7). In this context, WA^B = P ( X t r i a i = b\XCMRENT = a)x (2.9) P (Accept move? | X c u r r e n t = a & X t r i a i = h) The usual practice is to choose the trial moves symmetrically, so that P(XTV\3\ / ^ c u r r e n t = i) = P(XTRIA[ = i\XCURRENT = / ) . Then detailed balance can be satisfied by choosing the acceptance probabilities appropriately. One standard choice is accept a move from state i to state / with probability: minimum(l,p(/)/p(i)). Here, p(a) refers to the probability of finding the system in state a. An alternative choice for these probabilities is presented in section 2.2.4. Another important requirement, is that the probabilities be chosen so that the Markov Chain can access all states in S from any starting state. If these two conditions hold, then for sufficiently large N, XN should be essentially independent of the starting configuration xn.q Thus, recording the state of the Markov Chain after every N steps is equivalent to sampling the random variable XM repeatedly 3 . This is what is usually done in practice— sample a single Markov Chain every N steps. Continuous random variables Up to this point, I have discussed Markov Chains in the context of discrete random variables. This has served to make the discussion and theory rela-tively simple. However, these ideas can be extended to continuous random variables in a natural way. This topic will be discussed briefly in the current section, but continuous time/indexing variables and pathological functions will not be addressed. 2 To be explicit, P(Accept move?) is shorthand for P(Accept move? = yes) where "Accept move?" is a random variable with state state space {yes,no}. Similarly, P(Reject move?) is equivalent to P(Accept move? = no) 3 In fact, the desired condition is that the process is ergodic. For more on ergodicity see reference [30] 2. Computational Methods 11 Suppose we wish to construct a random variable X : f2 —> R, with probability density p{x). As before, the critical condition for generating a Markov Chain which converges to X is detailed balance. However in this context, the relevant states (subsets of the state space) are intervals (1-d volumes). That is, the condition of detailed balance can be written as: fbp(x)dx Jc p{x)dx In practice, one can usually carry out the procedure described above for discrete random variables by simply replacing the probability mass functions p with the probability density functions p. This corresponds to treating the state X = a G R, as the infinitesimal interval [a, a+dx]. Care must be taken, however, when using trial moves which map between unequal volumes of the state space. One common example of this is a scaling move, in which the current state is dilated or contracted. For the sake of concreteness, suppose that the initial state is i = x, and the trial move / is scaling by g. Then for detailed balance to be satisfied, we must have: ra+e rga+ge / p{x)dxWi^f= / p(x)dxWf^i (2.11) J o J ga ' which implies that p(a)W^f = gp(ga)Wf^ (2.12) Notice that one must include a factor of g when calculating the probability of accepting or rejecting this move. An application of this technique is discussed on page 124 of reference [31]. 2.2 Implementations of the Monte Carlo method and Applications to Physics In this section, I present two implementations of the M C technique for in-tegrating functions via importance sampling. As discussed in the previous section, the two main issues that must be decided are: i) how trial moves are to be determined, and ii) the probability of acceptance of the trial moves. Once these things have been determined, the algorithm outlined in section (B.4.2) becomes well-defined, and can be carried out mechanically. In order to make a good choice, one must consider the particulars of the integral in question. This discussion will focus on the ensemble average from statistical physics; specifically, as it pertains to the many-body problem. The 2. Computational Methods 12 implementations that will be presented are: i) standard Monte Carlo, and ii) Smart Monte Carlo (SMC). 2.2.1 Statistical physics and ensemble averages A central problem in statistical physics is the calculation of ensemble av-erages. Consider a system that can exist in many different states, and an observable, A, which is a function of the state. The ensemble average (A), is the value of A averaged over all configurations of the system, with each configuration properly weighted by its probability of occurrence. This can be written symbolically as: / A(w)dP(u>) (2.13) Jn where Q is the set of all possible configurations (the configuration space), and P is the corresponding probability measure function.4 It is often useful to identify the important degrees of freedom, {XQ, ..., Xjy} and integrate over these instead. Suppose that each Xi maps from the con-figuration space to the real numbers. Then the set X = {^ o> • • •, X^} cor-responds to a random vector with state space S C 1RN. For many physical systems, this random vector can be characterized by a probability density. If A can be written as a function of these degrees of freedom, then the ensemble average takes the standard form: i f A(x)p(x)dx (2.14) J s where p : S —> M is the (possibly unnormalized) probability density function, and Z is the normalization factor—a weighted sum of states, Z = f p(x)dx (2.15) Js that is called the partition function. The microcanonical and canonical ensembles will now be briefly dis-cussed. For simplicity, let us assume that X has been chosen so that its distribution is uniform if all states are equally likely.5 Let us also assume 4 See Appendix A , or references [32, 33] for more details. 5 This discussion is meant to be primarily illustrative, and will not address many of the subtleties of statistical mechanics. For a more complete account see, for example, references [34] or [35]. 2. Computational Methods 13 that the internal energy of a configuration, label by U, is a single-valued function of the state space S. Consider a system that is completely isolated, i.e. it does not exchange material, work or energy with the rest of the universe. If this "closed system" has a well-defined constant internal energy, then all states with this internal energy are equally likely. Symbolically, p(x) = <S(e(x) — U), where e(x) denotes the internal energy of the configuration corresponding to X = x. This is called the microcanonical ensemble, or distribution.6 Now suppose that the system is placed in a very large heat bath of constant temperature T, such that it may exchange energy with the rest of the universe, but nothing else. When it is in thermal equilibrium with the heat bath, the internal energy may fluctuate, but the average internal energy (U) will be constant. In this case, which is called the canonical ensemble, p(x) is given by the Boltzmann factor, p(x) = e _ / J e ( x ) (2.16) where /3 = 1/ksT, the reciprocal of the Boltzmann constant times the tem-perature. The Monte Carlo technique is well suited for the numerical calculation of many ensemble averages. The importance sampling technique is generally used, with a distribution function oip/Z for the sake of efficiency. 2.2.2 The many-body problem Monte Carlo simulations are frequently used to study systems of many par-ticles, such as liquids and gases. In addition, continuum models can often be discretized, and treated as collections of interacting particles. This is the technique that I used to model polymers. Suppose that the Hamiltonian of an iV-particle system can be decom-posed into kinetic and potential energy terms, with the potential energy dependent only on the positions. In this case, only the potential energy needs to be considered when calculating ensemble averages [36]. If in ad-dition, the observables of interest depend solely on particle positions, then these positions are the only relevant degrees of freedom. Let X be the ran-dom vector which represents positions of all of the particles, with state space S C M.Nd, where d is the dimension of the embedding space. Then as shown in the previous subsection, the ensemble average reduces to a standard Nd-dimensional integral over S. 0 However, not all closed systems have this d is t r ibut ion [34]. In some cases, certa in states of energy U may not be accessible f rom the system's init ial state. 2. Computational Methods 14 As always in the Monte Carlo approach, we integrate over the state space by constructing a Markov Chain. Each trial move usually consists of the displacement of a single particle. The particle to be moved can be chosen sequential (i.e. selecting each particle in turn), or randomly, with an equal probability of moving each particle. Alternatively, the trial move may consist of attempting to move all particles simultaneously. However, such a trial move is rarely used as the primary sampling technique for reasons of efficiency. Two common choices of trial moves and their corresponding acceptance probabilities will be presented in the next two subsections. 2.2.3 Standard Monte Carlo algorithm A standard trial move for moving a particle consists of choosing a new position with uniform probability within a "box" centered around the old position. The size and shape of this box may be selected to maximize code efficiency, or perhaps physical reasonableness. The shape is typically a cube or a sphere. The primary constraint is that the box boundary should have inversion symmetry, to ensure that given any two states, i and / , P(XTT[A\ = ^l^current — f) = -P(^trial = /[^current = This symmetry of moving between any two states allows us to choose a particularly simple acceptance probability. Suppose that the desired proba-bility density for the position random vector X is p(x). Then as mentioned in section B.4.2, a move from state i to state / , corresponding to a change in the particle positions from x ; to x/, 7 can be accepted with probability: PMC = P{Accept move? | ( A " c u r r e n t = i) & ( X t r i a i = /)) = minimum(l, p(x/)//9(xj)) (2-17) While the standard Monte Carlo (MC) method is fairly effective, there are situations in which other algorithmic settings are more appropriate. 2.2.4 Smart Monte Carlo algorithm The Smart Monte Carlo (SMC) algorithm is a simulation method that was developed by P. J. Rossky, J. D. Doll, and H. L. Friedman [36]. The basic idea is to use the equations of Brownian dynamics to determine the trial moves. By using information about the gradient of the potential8, this tech-nique is often able to produce convergence to the equilibrium state more 7 Or a change in positions infinitesimally close to these vectors 8 This is why the algorithm is "smart" 2. Computational Methods 15 quickly than the standard algorithm because of viscous damping. In addi-tion, the SMC method introduces an explicit time factor which facilitates the investigation of system dynamics. The motion of particles in a viscous medium is governed by the Langevin equation: r F + 5f r = -L + LJL£L (2.18) T TO where r is the position of the particle, 1/r is the frictional coefficient, TO is the particle mass, F is the deterministic force (the sum of the negative gradient of the potential due to all other particles plus any external force), and Sf is the random force. Note that we are ignoring long-range hydrodynamic interactions, and only considering the local (particle) drag. For further simplification, let us assume that the system is over-damped, so that r can be treated as zero, and also that r is a constant, independent of time and spatial position. Then integrating the equation of motion 2.18 yields the iterative step: Sr = j3FDSt + R (2.19) where Sr is the change in particle position in time step St, D is the diffusion constant for non-interacting particles, F is the mean value of the determin-istic force on the particle over this time interval, and R is the change in position due solely to the random force. Let R be a random variable with the Gaussian distribution: p(R) = {4DStn)-3/2exp{-R2/{4DSt)) (2.20) In addition, let us choose a value for DSt, and approximate F by the de-terministic force at the start of the time step. Then together with equation 2.19, we can define the trial move of a particle in the SMC simulation to be r i—> r + Sr. Now if DSt is sufficiently small, the motion determined by these trial moves can give a reasonable description of the time evolution of the system. However, even when St is large, Monte Carlo weighting can be utilized to ensure that the ensemble averages are correct. Consider a trial move that consists of moving a single particle from position Ti to rj. Let i and / be the initial and trial states of the system, with corresponding position vectors Xj and xy (recall that these vectors are instances of X, which describes the positions of all of the particles). Suppose that the desired probability density of the random vector X is p(x). Then 2. Computational Methods 16 in order to obtain detailed balance, we can use the acceptance probability: PSMC = mimmum(l, exp[(5r - PDStFi)2 - (-5r - (3D5tFf)2] x p(xi)/p(xf)) (2.21) where 5r = rj — r-j, and Fj is the force on the particle in state j. Notice that p(xi)/p(x.f) is simply the standard Monte Carlo weighting factor. The other part of this term is the "smart" factor. The SMC algorithm has been tested on a number of different models. It performed well, converging to equilibrium more quickly than the standard M C method in most situations [36]. For approximately the same amount of computational effort, the errors in the calculated averages were smaller with the SMC method than with the standard method. In addition, it has been demonstrated by Evans and Ritchie that the distribution of states under SMC evolution follows Smoluchowski theory as expected [37]. The extra computational overhead in the SMC method comes primarily from the calculation of forces. In many cases, the gradient of the potential has a simple form which is relatively inexpensive to calculate. However, if certain potential energy terms have "computationally nasty" gradients, the forces due to these terms can be ignored when carrying out the algorithm. In other words, while all potential energy terms must be included in the probability density function p(x), the gradients of unpleasant terms may be disregarded when determining the trial moves and when calculating the "smart" factor in the acceptance probability. Such a "hybrid" algorithm is acceptable as it preserves detailed balance. 2.3 Generation and Analysis of the Data In this section, the software and hardware which was used to generate and analyze data for the study is presented. Let us begin with the core of the project: the polymer simulation. 2.3.1 The Simulation Monte Carlo simulations were used to study the properties of a single poly-mer chain. The system was discretized and treated as a many-body problem, with each polymer represented as a collection of vertices. The random vector that determined the state space was constructed from the Cartesian coordi-nates of the vertex positions. Trial moves and acceptance probabilities were chosen according to a S M C / M C hybrid algorithm. A l l force terms were used, except for those related to an explicit bending stiffness. 2. Computational Methods 17 Each trial move consisted of moving a single vertex. Each "sweep" con-sisted of attempting to move each vertex sequentially. In some simulations, the polymer interacted with a movable plane. In these cases, the code at-tempted to move this plane once per sweep. Sampling was performed every M sweeps, where M was a large number which was determined from the simulation parameters. This number was chosen so that the sampled states were essentially statistically independent. Sampling consisted of making a number of "measurements" on the system, yielding specific instances of the random variables of interest. The first 10 samples of the system were not used in the calculation of averages, in order to give the system some time to reach equilibrium from its initial configuration. Generally, 200 samples were used to construct ensemble averages, with exceptions noted in the text. The simulation code was written in Digital Fortran, an extended version of the Fortran 90 language. The initial code was based on a SMC program that was developed in C++ by Ken Ritchie 9 and Evan Evans 1 0, which was in turn based on a M C program written by David Boa l 1 1 . Additional consul-tation for the code was provided by Carlos Camacho 1 2. A nearest-neighbor lookup table was used to increase the efficiency of the algorithm. Random numbers were generated using Numerical Recipes for Fortran 90 [38]. Specif-ically, the rani algorithm was used to calculate random numbers uniformly distributed from 0 to 1, and gasdev was used to generate random numbers with a Gaussian distribution. Eigenvalues of matrices were calculated with the IMSL library. The Digital Fortran Array Visualizer, along with its as-sociated libraries, was used for the three-dimensional visualization of the polymer. The simulation program was executed on Pentium II and Pentium III machines running Microsoft Windows. 2.3.2 Additional Analysis Much of the data analysis was performed by the simulation software during run-time. Additional analysis and graphing was carried out in Matlab, M i -crosoft Excel, and Microcal Origin. Errors in the ensemble averages were calculated statistically, using the standard error of the mean. This is sim-ply the standard deviation of the sample divided by the square root of the sample size. Non-linear fitting of the data was carried out using the 9 Department of Biological Science, University of Nagoya 1 0 Department of Biomedical Engineering, Boston University and the Department of Physics and Astronomy, University of British Columbia 1 1 Department of Physics, Simon Fraser University 1 2 Department of Biomedical Engineering, Boston University 2. Computational Methods 18 Levenberg-Marquardt method, as implemented in Origin. Further information about the models that were developed is presented in the following chapters. 3. Isolated Chain in Solution 19 3. ISOLATED CHAIN IN SOLUTION This chapter addresses the static properties of isolated, single chain poly-mers in solution. Some basic polymer theory is presented, with material drawn primarily from the text books by Doi and Edwards [18], Flory [39], and de Gennes [10]. The discussion focuses on three models: two standard analytical descriptions, and a computer representation which was designed to interpolate between the two. The primary objective of this chapter is to characterize the computer model, and to verify that it behaves as expected from theory. The chapter begins with a presentation of the different polymer models. This is followed by a study of size and shape properties in the large TV (or long polymer) limit. Concluding the chapter is an examination of the bending stiffness and the transition from stiff elastic rod to flexible chain behaviour. 3.1 Single Chain Polymer Models A polymer is a molecule which is composed of repeating units called monomers, The molecular composition of polymers can be very complex, with a wide variety of monomeric units, and multiple levels of structure. Although many different architectures are possible, this chapter focuses on the simplest case of an isolated, linear polymer in solution. I will study these systems using mesoscopic models, which do not include the detailed microscopic structure, but look at the larger scale features of the polymer chain. Two of the most common such models are the freely jointed chain (FJC), and the worm-like chain (WLC). 3.1.1 The freely join ted chain In the freely jointed chain model, the polymer is represented as a discrete chain of straight line segments. Furthermore, all of the line segments have the same fixed length, and there is no correlation between the directions of any of the line segments. In other words, the model is a random walk 3. Isolated Chain in Solution 20 with a fixed step/segment size. The restrictions on the length and angular distribution of the segments can be relaxed, to form more general random flight models. Another example of a random flight model is the "Diffusion Model" discussed in Chapter 4, Section 4.1.1. Formally, we can completely describe configurations of a F J C with N segments by a sequence of + 1 vectors, {R\, • • •, -Ryv+i}; where i?j repre-sents the position of the zth vertex. In some cases, the relative positions of the vertices is more useful, in which case we use the equivalent representa-tion {R, n , . . . , rjv}, where R = ( r i + . . . + rjy)/N is the centre of mass of the system, and r» = Ri+i — Ri is the relative position of vertex i + 1 with respect to vertex i. The F J C model is a good description of polymers that are highly flexible. We shall see in section 3.2 that most polymer chains with some flexibility can be accurately described by this model in the limit where they become very long (large N limit). Thus, the F J C model is convenient for calculating the universal properties of a large class of polymer models, due to its simplicity. 3.1.2 The worm-like chain The F J C model is not very good at characterizing polymers which have a non-negligible bending rigidity at the relevant length scales. One of the simplest models which captures the essence of such semi-flexible polymers is the worm-like chain. The worm-like chain model is sometimes called the Porod-Kratky model, as it was first presented by them in 1949 [40]. Essentially, it is a continuum model which includes an energetic penalty for bending (i.e. the polymer is represented by an elastic rod). The local bending energy per unit length is taken to be quadratic in the curvature. Formally, we represent the polymer as a smooth space curve T. Let r ( s ) be an arc length parameterization of T and let the contour length of the polymer be L. Then, the internal energy of the chain is given by the Hamiltonian: where kc characterizes the bending stiffness of the chain, and ^ is simply the local curvature1 of the chain. When kc is small, the model describes a 1 The basic geometric interpretation of the curvature in three-dimensional space is as follows. Consider a point, x, on the space curve, and the plane which best captures the space curve at this point. Now, consider the set of circles within this plane, whose 3. Isolated Chain in Solution 21 polymer which is very flexible. Alternatively, when kc is large, the worm-like chain is essentially a stiff rod. 2 In general, it is difficult to calculate the statistical features of the worm-like chain, as the path integral (or Green's function) corresponding to this energy functional does not have an easily solvable (e.g. Gaussian) form. However, there are a few useful quantities which can be calculated exactly, such as the end-to-end distance, and the persistence length. This will be dis-cussed further in subsequent sections. For additional information, including a summary of other known quantities, the review articles by E. Frey, K . Kroy, and J. Wilhelm are a good source [21] [41]. 3.1.3 Models with excluded volume and self-avoidance The two models above are examples of phantom chains—they are allowed to pass through themselves. In fact, there is no direct interaction between points which are sufficiently far apart on the chain. The lack of such "long-range" interactions makes these models easier to work with, but unrealistic for describing many physical systems. This problem can be remedied by introducing an excluded volume in-teraction to enforce self-avoidance. Kuhn introduced this approach in 1934 [42], and since then many different mathematical techniques have been de-veloped in order to study this interaction. While exact calculations on "self-avoiding" chains are quite difficult, a number of analytical results have been established. A summary of various approaches, including Flory's original estimation of the global polymer size [43], is presented by Doi and Edwards [18]. One particularly important result is that self-avoiding chains belong to a different universality class than phantom chains. However, if the interac-tion is "weak" (e.g. the excluded volume size is small relative to the step size), then this distinction may only be apparent for very large N (many polymer segments). Thus, a polymer with an excluded volume interaction may behave like a phantom chain at a given contour length, but if the length is increased sufficiently it will exhibit characteristically self-avoiding behaviour [22]. Physically, the excluded volume depends not only on the structure of perimeters pass through x. The osculating, or "kissing", circle, is the element of this set whose perimeter is the best approximation to the space curve at point x. The curvature is simply the reciprocal of the radius of the osculating circle. 2 Specifically, kc is small/large when the persistence length is small/large with respect to the contour length. This will be quantified more carefully in section 3.3. 3. Isolated Chain in Solution 22 the polymer, but also on the surrounding solution and its temperature. This effect is most pronounced in good solvents, and causes the polymer chain to swell. However, there is a temperature known as the G or Flory temperature, at which the repulsive steric effects are approximately balanced by attractive forces. At this temperature, the polymer looks almost ideal (i.e. as if it has no excluded volume). 3.1.4 The computer model Based upon the freely jointed chain and worm-like chain, I developed a com-puter model for use in Monte Carlo simulations. It was designed to approxi-mate and interpolate between these two analytical descriptions. Essentially, it was a random flight model with an explicit local bending stiffness that could be used to simulate a worm-like chain. The polymer was embedded in three-dimensional free space. The polymer was represented by a collection of spheres connected lin-early by perfect Hookean springs (they had harmonic potentials). These springs were in their lowest energy state when the centres of the spheres were separated by a specific bond length, bm. The spheres could be made "infinitely hard", which introduced an excluded volume interaction. An op-tional maximum bond length was included, which prevented the chain from passing through itself when an excluded volume was present. Finally, an energy penalty was added whenever there was a non-zero "exterior" angle between adjacent segments in the chain. (This is the angle which is sup-plementary to the vertex angle3. So if two adjacent segments fall upon the same line, i.e. there is no kink, the exterior angle is zero.) The model ne-glected intervertex (between the spheres) hydrodynamic interactions, and the attraction between spheres. This is a slight generalization of the Rouse model [44]. However, in the Rouse model, the bond length is always zero, and the excluded volume interaction is disregarded. In fact, the Smart Monte Carlo evolution of the chain was based upon Rouse dynamics, in which the hydrodynamic drag on each vertex is approximated by the drag on an isolated sphere. The sampling rate was based on the Rouse time. In other words, ac-cording to the pseudodynamical evolution of the system, the time between sampled configurations was given by the rotational relaxation time of the Rouse model: N2 h2 " 5 & <3-2> 3 The angle supplementary to 9 is 7r — 8 3. Isolated Chain in Solution 23 where Ny is the number of spheres, b is the root mean squared bond length, and D is the diffusion coefficient of each sphere. Preliminary simulations examined the relaxation time for a polymer initialized in a completely lin-ear state. Within a single Rouse time, the mean squared radius of gyration (which will be defined in subsequent sections) had decayed significantly be-yond the e _ 1 point, suggesting that the sampling rate was reasonable. How-ever, as the forces due to the bending stiffness were not used to generate trial moves, when the bending stiffness coefficient was non-zero the Rouse times and dynamic evolution were only approximate. The angular potential Recall that the worm-like chain has an internal energy term that depends on the line integral of the curvature squared. For the computer model, I approximated this by a finite sum, k Nv~2 ^bending = ^ ] fi (3-3) where Ny is the number of vertices (spheres) in the chain, and fi is a function of the configuration of the chain. Let r-j be the relative position of the i + 1th vertex with respect to the ith vertex, and Qi be the external angle of the i + 1th vertex. As the curvature is a local property, I let fi depend only on Oi, ri and ri+i. First, let us suppose that the polymer is basically inextensible (i.e. the springs between vertices are sufficiently stiff). Now from continuum me-chanics, we know that in the limit of small fluctuations (9 is sufficiently small), fi should be proportional to 2(1 — cos(0j)). Thus, a commonly used form for fi is the harmonic potential, fi ~ 92. This is quite good for small angles. However, when 8i approaches 7r, the curvature tends to infinity4, so we would like fi to approach infinity. The potential that I used with these two properties had the following angular dependence: / l ~ 4 t a n 2 | (3.4) Simple analysis shows that this angular bending potential has the desired form as 9 approaches 0, and as 6 approaches n. Now in my model, the bond lengths are not fixed, so fi must also have some length dependence. Let di be the average length of the vectors and 4 This is one interpretation of the model 3. Isolated Chain in Solution 24 r-j+i. That is, di = | r t l + 2 ' r t + l 1 (3.5) Since we are discretizing the integral in equation 3.1, we can write /j as K25S where «j represents the approximate curvature at vertex and Ssi is the approximate "length step" replacing ds in the integral. We can approximate the length step term by 5si = dj. Recall that the curvature is the reciprocal of the radius of the oscillating circle, and therefore must scale as 1/length.5 As it is a local property, I let KJ be proportional to 1/dj. Thus, fi was chosen to be proportional to 1/dj. So the final form of fi was simply the angular function given in equation 3.4 times the scale factor 1/dj. This yielded the angular potential: , ^ 2 t a n 2 | , , ^bending = 2fcc ^ ^ (3-6) i=l 1 Within the simulation, the trigonometric identity, tan 2 ^ = (1 — cos#)/(l + cos 9) was used to calculate the potential from the dot product of the relative position vectors. Units and scaled parameters To reduce the number of free parameters in the system, the physical param-eters were grouped together, and in some cases scaled to yield dimensionless "reduced" parameters. The unit of length was chosen to be the bond length (bm), with all positions and distances given as (dimensionless) multiples of this unit. The unit of energy was set to the thermal energy, UBT. Due to over-damping, the system did not have any mass dependence. Treating all units as unity, the spring constant for stretching the bonds, kt, and the bending stiffness coefficient, kc, were also dimensionless. Finally, the time step was multiplied by the diffusion constant D to give the single (dimen-sionless) parameter DT. (Explicitly, kt = kt • b^/ksT, kc = kclknTbm and DT = DT/b^, where kt, kc, and DT are dimensionful quantities.) Thus, the physically significant simulation parameters (and the symbols used to represent them) were: the number of vertices (including endpoints)(AV), the reduced radius of the spheres (Rs), the reduced intervertex spring coef-ficient (kt), the reduced bending stiffness coefficient (kc), and the reduced 5 Another way to think about it is that the curvature describes the change in angle per unit length, so in the discretized version it must scale as one over the finite length step. 3. Isolated Chain in Solution 25 step parameter (DT). Note that in most cases, instead of the number of vertices, I will talk of the number of segments, N = Ny — 1. Common parameter settings and final adjustments In the process of exploring the model, I settled on a few common parameter settings. First, the reduced step parameter was set to DT = 0.005. This produced results with reasonable efficiency. Second, two standard radius settings were chosen. To generate a chain which could not pass through itself, the maximum bond length had to be less than 2y/2 x Rs. Due to the "lurching" nature of Monte Carlo simulations, the maximum bond length was set to 2 \ /L9 X Rs. The hard-core sphere radius Rs was chosen to be 0.42045 bond length units. This made the mean of the maximum and minimum bond lengths approximately equal to one. The chain with this radius setting and the maximum bond length restriction enforced was called the "Self-Avoiding Chain", and will be henceforth referred to as SA. To generate a phantom chain, the long-range excluded volume inter-action was turned off. However, to make it locally consistent with SA, I enforced the minimum and maximum restrictions on the bond length con-sistent with a radius of 0.42045 (i.e. the bond length could not get smaller than 2 x 0.42045, or larger than 2 \ /L9 X 0.42045). This model was called the "Partial Phantom Chain", and labelled by P P C . In all other cases, the maximum bond length restriction was turned off. As for the intervertex spring constant, in most simulations it was set to kt = 100. This is the default, unless otherwise specified. To make the kt setting explicit, it was sometimes appended to the label P P C or SA (e.g. PPC100 for a partial phantom chain with kt = 100). 3.2 Shapes and dimensions in large N limit The static properties of the polymer models were studied in the large N limit. In this regime, there are some well-established universal properties which are independent of the local structure of the chain. The results of the computer simulation were compared with the expected behaviour. It was found that the computation model was quite reasonable. 3. Isolated Chain in Solution 26 3.2.1 Global length scales One common way of characterizing the mean global length scale of linear polymer is the root mean squared end-to-end distance, \J(r|e), where r e e is simply the distance between the endpoints. For the FJC, the mean squared end-to-end distance is easy to compute [18, 39]. For a FJC with N segments, JV N i=i j=i N i=\ i>j Since there is no correlation between r-j and Tj for all i ^ j , the second collection of terms is zero. Finally, since the bond length is constant, we have for all N: (r2ee) = Nb2 (3.8) If the bond length is not constant, but there are still no angular correlations, then b is simply the root mean squared bond length. Kuhn [45, 46], made a correspondence between this model and the large N limit of "real" polymers with no long-range interactions by the following argument. We can describe a polymer by a collection of N relative posi-tion vectors, {ri} and its centre of mass (i.e. as a generalized random flight model). Provided there are no long-range interactions and the polymer has some flexibility, the angular correlation between rj and ri+k must approach zero as k approaches infinity. In fact, ^2f=i+iiri ' rj) must converge, and this sum should be independent of i, provided that segment i is not near to the endpoints. Then we see that {r2e) in equation 3.7 is proportional to N, since the first set of terms is simply Nb2ms, where bTms is the root mean squared bond length. Thus, we can make an equivalent F J C by choosing the Kuhn statistical segment length bpc, and the statistical number of segments NK such that: bkNk = L (3.9) where L is the contour length, or maximum end-to-end distance, and (3.10) 3. Isolated Chain in Solution 27 It should be noted that the term "equivalent bond length" is sometimes used to refer to ^/{r2e)/N. For the worm-like chain, it has been shown that the dimensionless Kuhn length is twice the dimensionless bending stiffness (kc) in the limit of long chains. I will discuss this further in section 3.3. This large N correspondence with the F J C chain is not restricted to (r 2 e ) . In fact, all static, global shape properties of polymers with no long-range in-teractions can be described by a F J C chain, with the local properties simply serving to establish the equivalent bond length [18]. Doi and Edwards [18] also provide a way to think about this correspondence. Consider dividing the polymer into M = N/k submolecules, each consisting of k segments. Provided k is sufficiently large, the directions of these submolecules will be uncorrelated, and the polymer will act as a freely jointed chain, provided M is much larger than k. Another measurement of the global length scale is the radius of gyration, rg. It is simply the root mean squared distance of the polymer segments (or in the discrete case, of the vertices) from the centre of mass. As calculated in [39] for a FJC , 2 livfr 2(/V + 2) {ra'~Q N+l [6AS-} Hence, as shown by Debye [47], in the large N limit (r2g) = <d>/6 (3.12) Again, this limiting case generically describes long, single chain polymers with no long-range interactions. For the worm-like chain, the correspondence can be shown explicitly [48]. While it is difficult to calculate anything exactly in the presence of an excluded volume interaction, the scaling law for (r 2 e) has been estimated by various theoretical approaches. The generally accepted value for the scaling law in three dimensions is (r2ee) ~ TV1-2 (3.13) as derived from an argument by Flory [43]. In fact, this results may be generalized to other dimensionalities as: <r e 2 e )~/V 2 ", ^ = 3/(^ + 2) (3.14) where d is the number of dimensions of the embedding space [49]. This formula is very accurate for up to four dimensions. Above four dimensions, the chain basically looks ideal (like a phantom chain) with v = 1/2. For more details of this derivation, in addition to some other approaches see [10, 18]. 3. Isolated Chain in Solution 28 Results and analysis of computer simulations In order to verify that the computer model worked as I expected, I ran a number of simulations on the contour length, mean squared end-to-end distance, and mean squared radius of gyration, as a function of the number of segments TV. Furthermore, these simulations were performed for a number of different parameter settings. Due to computational limitations, TV was only increased to 90 segments. First, I examined the two principal parameter sets: the partial phantom chain P P C and the self avoiding chain SA. The intervertex spring constant was set to kt = 100, in order to make the bonds stiff enough to resemble a freely jointed chain. The bending stiffness coefficient kc was set to zero. Various properties were examined as a function of TV. As shown in Figure 3.1, the contour length grew linearly with the number of segments, as expected. The results of a weighted linear regression to the data are also presented in Table 3.1 (data was weighted by 1/error2). Due to the variable length of the intervertex bonds, the constant of proportionality was slightly larger than 1 (by about 1%), but this was acceptable. The vertical offset was negligible. Model A (Offset) B (Prefactor) PPC100, kc = 0 0.005(5) 1.0116(2) SA100 kc = 0 -0.008(4) 1.0141(2) PPCO kc = 0 -0.004(5) 1.0168(2) SAO kc = 0 -0.015(5) 1.0206(3) Tab. 3.1: Contour length (L) fit for various models: L = A + B x TV. Partial Phantom Chains (PPC) and Self-avoiding Chain (SA) with reduced tether stiffness values of kt = 100 and kt = 0 were examined. There was no explicit bending stiffness (kc = 0). Figure 3.2 shows (r 2 e) and (r 2) as a function of the number of segments for these two chains. The graphs were fairly linear on the log-log plot, for TV approximately greater than 10, indicating power law scaling for large TV. A power law fit was performed on the data for TV > TVi , where TVi was chosen to be within the apparent large TV regime (region of linearity in the plot). The results are presented on the graphs, and summarized in Table 3.2 6 . For the phantom chain, (r 2 e) and (r 2) both scaled approximately linearly with TV as expected. For the self-avoiding chain, (r-2e) scaled as T V 1 - 2 3 ' 2 ) and (r 2) scaled as T V 1 - 2 3 ' 1 ) , which was slightly larger, but comparable to the expected 6 The number in brackets gives the estimated uncertainty in the least significant digit. 3. Isolated Chain in Solution 29 Number of segments Number of segments Fig. 3.1: Mean contour length vs. number of segments for the computer model. (a) Partial Phantom Chain with reduced tether stiffness of kt = 100 (PPC100), (b) Self-avoiding Chain with reduced tether stiffness of kt = 100 (SA100). In this simulation, there was no explicit bending stiffness, kc = 0. The results of a linear fit are superimposed. 3. Isolated Chain in Solution 30 scaling of i V 1 - 2 . Furthermore, the prefactors were approximately as expected for both of the models: close to unity for (r2e) and around 1/6 for (r 2). Model (•> .2 v eel (r2g) Ni a b a b F J C 1 1 [39] 1/6 = 0.167 [47] 1 Flory (SA) 1.2 [43] PPC100, kc = 0 0.92(4) 1.03(2) 0.194(6) 0.97(1) 6 SA100, kc = 0 1.09(6) 1.23(2) 0.176(6) 1.23(1) 10 PPCO, kc = 0 1.02(6) 1.00(2) 0.212(8) 0.95(1) 6 SAO, kc = 0 1.1(1) 1.22(3) 0.18(1) 1.23(2) 13 Tab. 3.2: Scaling laws for (r^e) and (r^) for various models. Da ta was fit to y = axb for N between N1 and 90 (inclusive). Parameters used in the PPCfOO, SA100, P P C O , SAO simulations are noted in the legend of Figure 3.1. Intervertex spring coefficient To study the effect of the intervertex spring coefficient, kt was set to zero, and the same measurements made. The results were very similar to the chains with kt = 100, as can be seen in Tables 3.1 and 3.2. First, the mean contour lengths still grew linearly with N, though the coefficient was about 2% greater than unity, due to variations in the bond lengths. However, this difference is not drastic. For the P P C chain the scaling laws differed from the stiffer linking potential by ~3 %, and the prefactors by about ~10%. For the SA chain, the differences were smaller. For the mean squared radius of gyration, the values of the power law fit were identical, within error. For the mean squared end-to-end distance, the scaling law exponent differed by less than 1%, and the prefactors by less than 5%. So in summary, with the maximum and minimum cut-offs in place, the magnitude of kt did not seem to have a gross effect on the average global size of the polymer. In the absence of a bond length cut-off, however, preliminary tests indicated that the polymer became quite stretched out, on average. Excluded volume Next, to investigate the excluded volume effect, simula-tions were run for spheres with the following hard-core sphere radii: 0.005, 0.1, 0.2, 0.25, 0.3, 0.35, and 0.4. kt was set to 100. The maximum bond length was not turned on, as 1 did not want to enforce strict self-avoidance. (Although for sphere radii greater than l/(2\/2), exploratory simulations indicated that the maximum bond cut-off did not have a large effect when 3. Isolated Chain in Solution 31 ^ •— •— • i •— •— • 1 10 100 1 10 100 Number of segments Number of segments Fig. 3.2: (r1e) and (r^) vs. number of segments for various models. The data is presented in a log-log plot. A fit to the equation y = axb was performed for N > Ni, and the results superimposed. (a,c): P P C 1 0 0 , kc = 0, Ni = 6. (b,d): SA100 kc = 0,Ni= 10. 3. Isolated Chain in Solution 0.2 0.3 0.4 Sphere radius 0.2 0.3 0.4 Sphere radius Fig. 3.3: Mean squared end-to-end distance prefactors and exponents vs. sphere radius. Various fits of ree2 = aNb for the parameters a and b were per-formed for various sphere size. In these simulations, kt = 100, and kc = 0. The results are summarized here. the linking stiffness was this large.) The power law scaling of (r 2 e) with the number of segments was investigated for each sphere size, for N up to 89. A summary of the scaling law exponents and prefactors, including the results from the P P C and SA chains, is presented in Figure 3.3. The prefactors were approximately unity, which suggested that the sphere size did not significantly change the effective bond length within the range of investigation. However, the scaling exponents varied from around 1 for the FJC , to slightly more than 1.2 for the largest sphere sizes. As men-tioned in subsection 3.1.3, the continuum ball and spring model should also exhibit a sharp distinction between phantom and excluded volume classes for sufficiently large A'' (see reference [22]). Thus, A'' = 89 clearly did not represent the large N limit for any of the sphere sizes, save the largest two (which fortunately included the SA chain at a sphere size of 0.42045). Explicit bending stiffness Next, I studied the large N behaviour of the model for different values of the bending stiffness coefficient kc. I calculated (r2e) and (r2) as a function of the number of segments for numerous values 3. Isolated Chain in Solution 33 of kc from 1 to 9. The intervertex stiffness was held at kt = 100 First, I examined the behaviour of the P P C model by plotting the values on a linear scale. The graphs were clearly concave up (had a positive second derivative) for small N, then appeared to approach a linear regime. However, within the range examined, only the graphs for kc € {1,2} had a linear region that was long enough to attempt a determination of the power law scaling. In fact, for the largest values of kc, it was not obvious whether we had entered the large N regime. As presented in Figure 3.4, (r2e) and (r2) were plotted on a log-log scale, and a power law fit was performed for N > Nj. The results are summarized in Table 3.3. The scaling coefficients were all roughly unity within error, showing the FJC-like behaviour of the model. The prefactors were larger than the expected values of 2 x kc, suggesting that the "effective" bending stiffness may have been larger than expected. Furthermore, the prefactors for the kc = 2 model were not quite twice as large as for kc = 1: approximately 70 % larger for (r 2 e) and 50% larger for (r2). However, the errors in the prefactors were so large, that this was not a concern. Since N was not large enough to fit the data over even a single decade, the numerical results were not persuasive. Indeed, this test was not intended to establish the scaling law, but as a quick check of the reasonableness of the model. Model (rle) Wi a b a b PPC100, kc = 1 3.5(7) 1.04(6) 0.46(6) 1.08(4) 20 PPC100, kc = 2 6(2) 1.00(9) 0.7(2) 1.07(8) 30 SA100, kc = 1 5(2) 0.98(8) 0.51(9) 1.10(4) 30 Tab. 3.3: Scaling laws for (r^e) and (r%) for chains wi th an explicit bending stiffness. Da ta was fit to y = axb for N between Ni and 90 (inclusive). In order to check the reasonableness of the prefactors for other values of kc, linear fits of (r 2 e) and (r2) versus N were performed on all of the P P C data, for N > 40. The fit was constrained to pass through the origin. A sample of the data is presented in Figure 3.5. The results are summarized in Figure 3.6. The prefactors for kc ^ 0 were generally larger than the predicted values of kc x 2 for (r 2 e), and kc/3 for (r2). For values kc from 1 to 4, the prefactor for (r2e) grew roughly linearly with a proportionality constant of 2 and an offset of 1. This suggested that the "effective" bending coefficient may have been slightly larger than the value of kc put into the model. Within this range of kc, the prefactors for (r 2) were slightly less 3. Isolated Chain in Solution 34 Number of segments 10 100 Number of segments Fig. 3.4: (r1e) and (r2,) as a function of segment number (N) for two part ial phan-tom chains ( P P C ) wi th non-zero explicit bending stiffness. The data is presented in a log-log plot. A fit to the equation y = axb was performed for N > Ni, and the results superimposed. (b,d) P P C 1 0 0 , kc = 2, Ni = 30. (a,c): P P C 1 0 0 , kc = 1, N\ = 20. 3. Isolated Chain in Solution 35 Number of segments Number of segments Fig. 3.5: (ri) and (r2q) for P P C 1 0 0 chain wi th kc = 6. A linear fit was performed on this data for ./V > 40, and the results superimposed. than one sixth of the ( r 2 e ) prefactors. For larger values of kc, the prefactors of both quantities did not increase as quickly as predicted by the analytical model, "rolling" below linearity. This indicated that the range of A'' analyzed was not yet in the large N regime for these values of kc. Next, I modelled the S A chain. Only for kc = 1 was there a large enough linear region in the log-log plots of ( r 2 e ) and ( r 2 ) to attempt a power law fit. However, the region was still too small to give more than a vague indication of the relevant scaling laws. The graphs are presented in Figure 3.7, and the fitting values displayed in Table 3.3. The scaling exponent for (ree) w a s nearly unity within error, rather than the characteristic value of approximately 1.2 for self-avoiding chains. One interpretation is this: for an isolated worm-like chain with a finite radial thickness (excluded volume), increasing the length and the bending stiffness is equivalent to decreasing the radial thickness, as there are only three length scales. Thus, increasing kc effectively decreased the excluded volume effect. To investigate this further, more simulations could be performed for larger values of A^. In section 3.3,1 will investigate the effective bending stiffness and the reasonableness of the model by looking at angular correlations. 3. Isolated Chain in Solution 36 O CD CD i_ Q. m 20 18 16 14 12 10 8 6 4 2 0 • <r 2> ee A <r 2> * 6 g " A : I 5 2 4 6 8 Bending stiffness coefficient 10 Fig. 3.6: Prefactor of linear fit to ( r^ e ) vs. N for various settings of the bending stiffness kc. The simulations were of a P P C 100 chain. 3. Isolated Chain in Solution 37 Fig. 3.7: (rle) and (r2g) as a function of segment number (./V) for SA100 chain with kc = 1. The data is presented in a log-log plot. A fit to the equation y = axb was performed for N > 20, and the results superimposed. 3. Isolated Chain in Solution 38 3.2.2 Shape properties In addition to global length scales, I investigated the average shape of poly-mers in the large N limit by studying the quantities: the asphericity, S3, and the principal moments of inertia. The asphericity and S 3 are based on the radius of gyration tensor [50]: Qa/3 = 53 [Xi<a ~ xiA\xi,P ~ XjA (3.15) where Ny is the number of particles/vertices, a,(3 G {1,2,... ,d} where d is the number of dimension in the embedding space, and Xi,a represents the ath component of the position vector for the ith vertex. For example, (Xi,i,X\t2,X\,$) represents the position of the first vertex in three dimen-sional space. Qap has d eigenvalues, which shall be labelled by qa, where a is an index that runs from 1 to d as before. Let q be the mean eigenvalue. Then, the familiar radius of gyration is simply dq. The asphericity measures how much the shape deviates from a sphere. It is defined as [51, 52]: -^ssnnjEV (316) v ' Q=l It varies in value from 0 for a spherical distribution, to 1 for a completely collinear distribution. The quantity S 3 allows one to distinguish between prolate (football-like) and oblate (pancake-like) shapes in three dimensions. It is defined as [52]: S3 = ( 3 . 1 7 ) q6 Positive values of S 3 indicates a prolate shape; negative values of S 3 indicate an oblate shape. For more information about these quantities, see reference [53]. Another standard way to study the mass distribution of an object is by examining the principal moments of inertia. These (scalar) moments are simply the eigenvalues of the moment of inertia tensor [54]: NV lap = J2 m(Xl~f6*P - Xi,«Xi,n) (3-18) i=l 3. Isolated Chain in Solution 39 where 5an is the unit tensor (6ap = 0 if a ^ (3, 5aa = 1), summation is implied over repeated indices (X? = 2 ^ 7 = 1 ^ 7 ) 1 a n d the rest of the symbols are as defined in equation 3.15. In three dimensions, there are three principle moments of inertia. In order from smallest to greatest, let them be labelled by I\, I2, and I3. As I was interested in the average shape, I divided the two smaller principle moments of inertia by the third, for each configuration sampled. These two dimensionless ratios, I1/I3 and /2/T3, were then averaged over all configurations. Information on the magnitude of the principle moments of inertia is given by the radius of gyration. I expected the average shape of the polymer to be independent of N, provided that N was sufficiently large. In the absence of long-range inter-actions, the asphericity, and S 3 , have already been shown to approach such universal values [53, 55]. In fact, these quantities have been calculated ana-lytically [53, 55]. An estimation of (.A3) in the presence of excluded volume has also been obtained via an e = 4 — d expansion [53]. These analytical results are displayed in Table 3.4 for comparison with the simulation data. I am not aware of any published analytical calculations of the mean ratio of principle moments of inertia. Model (S3) (h/h) (h/h) Analytical phantom chain PPC100, kc = 0 PPC100, kc = 1 PPC100, kc = 2 0.394 [55] 0.385(3) 0.41(1) 0.43(1) 0.475 [53] 0.453(7) 0.49(3) 0.53(2) 0.306(3) 0.290(8) 0.278(9) 0.868(2) 0.863(5) 0.863(6) Analytical SA chain SA100, kc = 0 SA100, kc = 1 0.416 [53] 0.451(3) 0.43(1) 0.575(8) 0.54(3) 0.260(2) 0.27(1) 0.871(2) 0.867(5) Rs = 0.2, kc = 0 0.415(4) 0.512(8) 0.286(3) 0.871(1) Tab. 3.4: Large N shape properties for various models. The properties presented are the mean asphericity ((.A3)), {S3), and the mean ratios of the three principal moments of inertia ((I1/I3) and (J2//3)) Results and analysis of computer simulations Shape properties were calculated for these models with N ranging from 3 to 90. The P P C and SA chains were studied for various settings of the bending stiffness coefficient. Additionally, a polymer with a hardcore radius of Rs = 0.2 and no explicit bending stiffness, was simulated to investigate the dependence on excluded volume. A l l of these simulations had a intervertex 3. Isolated Chain in Solution 40 (b) S A 1 0 0 20 40 60 80 N u m b e r of s e g m e n t s I N u m b e r of s e g m e n t s Fig. 3.8: Mean asphericity vs. number of segments for various models, (a) P P C 1 0 0 , kc = 0, (b) SA100 kc = 0. fn order to determine the approximate large N l imi t , a linear fit w i th slope 0 was performed on each data set from N = 20 to N = 90. The results are overlaid. spring constant of kt = 100. For all of the models with no explicit bending stiffness, the shape prop-erties appeared to approach stable, constant values within the range of N investigated. Numerical values for the large N shape properties were deter-mined by performing a weighted linear regression on this data for N > 20, with the slope constrained to be zero. An example of the data and fitting is presented in Figures 3.8, 3.9, and 3.10. The results are summarized in Table 3.4. When kc was non-zero, only a few of the models seemed to enter the large ./V regime. As with (r 2 e) and (r 2), the only models worth measuring were kc £ {1,2} for P P C , and with kc = 1 for SA. Averaging was performed on this data for N > 70, to get a rough estimate of the large N values. The results of the fit are included with the others in Table 3.4. The average shape of all of the models was prolate, consistent with the results of previous studies. This was implied by (S3) > 0 for all the models. In addition, (h/h) always ranged from 0.86 to 0.88, while (I\/h) generally took values around 0.3. This indicated that the two largest principle mo-3. Isolated Chain in Solution 41 Fig. 3.9: Mean S3 vs. number of segments for various models, (a) PPC100, kc = 0, (b) SA100 kc = 0. In order to determine the approximate large N limit, a linear fit with slope 0 was performed on each data set from = 20 to N = 90. The results are overlaid. 3. Isolated Chain in Solution 42 (a) P P C 1 0 0 —i— 60 - 1— 80 Number of segments co 100 (b) S A 1 0 0 - 1— 20 — I — 40 Number of segments Fig. 3.10: Mean ratios of principle moments of inertia vs. number of segments for various models, (a) P P C 1 0 0 , kc = 0, (b) SA100 kc = 0. For each configuration, the three principle moments of inertia were calculated, and ordered from smallest to largest as I\, I2 and I3. Circular points represent (I1/I3); square points represent (12/13)- In order to determine the approximate large N l imit , a linear fit wi th slope 0 was performed on each data set from N = 20 to N = 90. The results are overlaid. 3. Isolated Chain in Solution 43 merits of inertia were usually fairly close in magnitude, while the third was usually much smaller. This is consistent with a prolate shape. In addition, it suggests that the average shape had roughly cylindrical symmetry, though more measurements would be required to quantify this. The results for the PPC100 phantom chain with kc = 0 were quite close to the analytical prediction (less than 3% smaller). As kc was increased but the contour length held fixed, the polymer became increasingly prolate and asymmetrical. For sufficiently large N it was expected that the shape would be independent of kc. This suggests that N was not large enough to put the system into the large N regime for these values of the bending stiffness. The excluded volume interaction also caused the polymer to become more prolate and asymmetrical. This effect was greater than increasing the bending stiffness to kc = 2 in the PPC100 chain. The mean asphericity was approximately 8% larger than the analytical predictions. Again, this may have been due to insufficiently large N. When the bending stiffness was in-creased in the self-avoiding chain, the mean values of A% and S 3 decreased. This may have been due to the corresponding decrease in the effective ex-cluded volume. For the model with an excluded volume given by Rs = 0.2, the val-ues of each shape property fell in-between the corresponding P P C and SA values. Thus, within this range of N, the model with a sphere radius of 0.2 did not show universal self-avoiding or phantom chain behaviour. This was consistent with the results for the SA100 chain with non-zero bending stiffness. 3.3 Persistence length and cross-over from elastic rod to freely jointed chain behaviour For a polymer with some bending rigidity, as its length increases, there is a transition from a rod-like structure to a randomized structure which looks like a freely jointed chain. In fact, this behaviour can be seen in many of the figures in the previous section. For example, in Figure 3.5, the graphs of (r-gg) and (r 2) look roughly quadratic for small N, before entering the linear regime. I will examine this cross-over for the worm-like chain, and compare the behaviour to the results of the computer simulation. Consider the tangent-tangent correlation function7 of the worm-like chain. Due to thermal fluctuation, it decays exponentially with length [34]: (t(s) -t(s')) = e-'8-8''/'" (3.19) 7 Alternatively, the angular correlation function 3. Isolated Chain in Solution 44 where t(s) is an arc-length parameterization of the unit tangent vector. The decay constant lp is called the persistence length. It is an important measure of the stiffness of the polymer, and defines an intrinsic length scale.8 In three-dimensions, lp has been calculated exactly as [34, 39, 40, 56]: lP = kc/(kBT) (3.20) From this result, another important quantity can be calculated: the mean-squared end-to-end length, (r2e). This calculation appears in a number of places, [34, 39, 40, 56], yielding: (r2ee)=2l2p(e-L^-l + L/lp) (3.21) where as before, L is the contour length of the worm-like chain. Two limiting cases are worth noting. First, when the polymer becomes very stiff, i.e. L <C lp , (r2e) approaches L2{1 — L/(3/ p )) , and ultimately I?. So the polymer looks like a rigid rod, as expected. Secondly, when the polymer is very flexible, i.e. lp <C L, (r2e) approaches 2lpL = (2lp)2 (L / (2lp)), so the polymer looks like a random walk, or freely jointed chain, with Kuhn length 2lp. The mean squared radius of gyration has also been calculated [48]: (r2g)=l2p((r2ee)/L2-l+L/(3lp)) (3.22) As expected, when L becomes much longer than the persistence length, (r2) approaches (l/6)(2lp)2(L/(2lp)), consistent with a F J C with Kuhn length 2lp. Note that the transition of both {r2e) and (r2) to F J C behaviour is characterized by an exponential decay with decay constant lp. In the next section, I quantify the cross-over behaviour of the computer model by cal-culating lp as a function of the kc simulation parameter. 3.3.1 Results and Analysis of Computer Simulations For this set of simulations, I specialized to the P P C and SA parameter sets, with an intervertex spring constant of kt = 100. The bending stiffness 8 A n alternative definition of the persistence length is given in references [39, 40, 56] for polymer models which consist of a finite number of segments (general random flight models). Consider some segment i far from the ends of the chain. The persistence length may be defined as the average sum of the projections of all subsequent segments (j > i) on segment i, in the limit that the chain becomes infinitely long. For the worm-like chain, this definition agrees with the one given in the text [39, 40, 56]. In the present study, I shall treat and measure lp as a decay constant. 3. Isolated Chain in Solution 45 coefficient kc was varied from 1 to 10 in integer steps. The tangent-tangent correlation function was then calculated as a function of the number of segments separating the tangents. As kt was fairly large, the "separation number" was approximately equal to the arc-length separation distance in our system of units. As expected, when kc was non-zero, (t(i) • t(j)} exhibited roughly expo-nential decay with the separation number \i — j\. The data was fitted with a single exponential, e Jt/*i ^ and the decay constant t\ was identified as the persistence length, lp. This in turn, defined an "effective bending stiffness coefficient" fcefj = lp (within our system of units), according to equation 3.20 for the worm-like chain. A sample of the data is presented in Figure 3.11. For the P P C model with kc — 0, there were no apparent angular cor-relations between the directions of the segments. In this regard, it truly looked like a freely jointed chain. However, for the corresponding SA model with kc = 0, short-range angular correlations were apparent. A fit to the exponential yielded a decay length of ti = 0.82(3). This effect was not sur-prising, as the excluded volume prevented the polymer from doubling back on itself. At non-zero temperatures, the steric barrier becomes "smeared out", resulting in an effective bending rigidity. The persistence length as a function of kc is presented in Figure 3.12 for P P C , and Figure 3.13 for SA. As expected, these graphs were fairly linear. A linear regression was performed on the data, y = A + B x x. This allowed the effective bending stiffness coefficient, fceff, to be determined as a function of the simulation parameter kc. For P P C , the offset, A , was 0.85(9), and the slope, B, was 1.09(4). So the angular potential that was used caused the persistence length to grow slightly more quickly than expected. In addition, it did not extrapolate perfectly to kc = 0. Overall, though, the model was reasonable, and the deviation of the angular correlation function from ideal worm-like chain behaviour was at least fairly systematic. In addition, the results were consistent with the (r2e) prefactor data in Figure 3.6, which indicated that fceff grew approximately linearly with kc, with an offset of about 1. For SA, the offset was given by A = 1.0(1), and the slope was fitted to be B = 1.12(4). It was not surprising that the slope was slightly higher in the presence of an excluded volume interaction. Furthermore, the non-zero bending stiffness for kc = 0 was consistent with expectations. For the largest values of kc, moderate fluctuations arose in the angular correlation data as separations between segments became large (see Figure 3.11 (d)). In particular, small wiggles followed the rapid decay at short sep-arations. Fortunately, the measured decay constants were quite reasonable within the range of kc tested despite these fluctuations. To investigate this 3. Isolated Chain in Solution 46 Fig. 3.11: Tangent-tangent correlation function ((t(i) -t(j))) vs. separation number (|i - j\) for. various chains, (a) PPCfOO, kc = 0, (b) SA100, kc = 0, (c) P P C 1 0 0 , kc = 1, (d) PPCfOO, kc = 9. Cha in length was 89 seg-ments; tethering stiffness was kt = f 00. The data was fitted to a single exponential, y = ex/n, and the results superimposed above. 3. Isolated Chain in Solution 47 Fig. 3.12: Persistence length as a function of the bending stiffness kc for PPC100. Persistence length was determined from the decay constant (ti) in the fitting equation (t(i)4(j)) = e _ ' l _ J ' / t l for the tangent-tangent correlation function. The results of a weighted linear regression are overlaid. 3. Isolated Chain in Solution 48 L inea r R e g r e s s i o n : 1 ' 1 ' I ' I ' I ' 1 ' 1 ' 1 ' 1 ' 1 1 I ' 1 0 1 2 3 4 5 6 7 8 9 10 11 Bending stiffness coefficient Fig. 3.13: Persistence length as a function of the bending stiffness kc for SA100. Persistence length was determined from the decay length (ti) in the fit-ting equation (t(i) • t(j)) = e~l l - J ' /* 1 for the tangent-tangent correlation function. The results of a weighted linear regression are overlaid. 3. Isolated Chain in Solution 49 further, simulations were performed for much larger kc. It was found that linearity in the decay constants was lost between kc = 15 and kc = 20. One possible explanation is that the time between samples was not long enough in this extreme stiffness regime. Since angular forces were not explicitly calculated for trial moves, the time evolution of the system did not follow Smoluchowski theory when these forces dominated. Thus, it would not be surprising if the estimated Rouse time became a poor indicator of the time between independent configurations in this regime. Nevertheless, the overall behaviour of the computer model was acceptable for both large and small values of N, and for kc from 0 to 10. With the appropriate parameter set-tings, the model could be made to behave like a worm-like chain or a freely jointed chain, as desired. 4. Single Chain Tethered to an Interface 50 4. S I N G L E C H A I N T E T H E R E D T O A N I N T E R F A C E Having established that the computer model behaved as expected for an iso-lated chain in solution, I now describe results for a slightly more complicated system: a single chain tethered at each end to an interface. This system is known as a polymer loop. The chapter is divided into three sections. First, two models for studying polymer loops are presented, one computational, the other analytical. Next, the shape properties of polymer loops as a func-tion of contour length are investigated, primarily in the "large N" (i.e. long polymer) limit. Scaling laws are derived from the computational results, and these relationships are compared to the analytical theory. The chap-ter concludes with an examination of the effect of bending stiffness on the shape of a polymer loop. In this section, the computational worm-like chain is compared not only to the analytical model, but also to the equivalent freely jointed chain. 4.1 Single Polymer Loop Models The specific system investigated in this chapter is a single linear polymer tethered at each end to an infinite impenetrable plane. In other words, the polymer is restricted to one side of the plane, with each end constrained to remain a fixed distance from the plane. No additional constraints are im-posed on the ends, so they are free to "skate" along the surface, which implies a stress-free lateral boundary condition on the loop. The zero-temperature interaction potential between the polymer and the surface is simply an in-finite square barrier, so adsorption is not considered. Finally, the ends are tethered as close to the plane as possible. (Clarification of this statement will be provided in the following subsections.) This system has been studied in the literature, most notably in [7] and [17]. In the following subsection, I present an analytical model which is discussed in these papers. 4. Single Chain Tethered to an Interface 51 4.1.1 The Diffusion Model A useful analytical description of a flexible chain molecule is given by the diffusion equation. In particular, a polymer is modeled as the path taken by a particle as it undergoes Brownian motion. As mentioned in Chapter 3, this is an example of a random flight model. Although the model is somewhat limited, as it only describes flexible polymers1, it is very useful, due to the ease with which its properties can be calculated. I shall refer to it as the Diffusion Model (DM). Another name for this model is the Gaussian chain [18], due to the form of its distribution function, which shall be discussed shortly. It has also been called the random flight model. The differential equation which describes the probability density of the Diffusion Model is: a G ( r9r ; 5 ' 5 / ) - \ V 2 G ( r , r', s, s>) = 5(r - r')5(s - s') (4.1) This is simply the standard diffusion equation (also known as the heat con-duction equation), with the diffusion coefficient set to 6/6, the time pa-rameter replaced by s, and a delta-function source (i.e, the "particle" is localized at r' when s = s'). To make the correspondence with a random flight polymer model complete, b is interpreted as the segment length, and the parameters s and s' are said to represent distance as measured along the polymer's contour. Then G(r,r',s,s') represents the probability density of finding "segment s" 2 at position r, given that "segment s'" is at position r'. (That is, G(r, r', s, s')dV is the probability of finding the segment within the volume element dV centered about r.) Another way of putting it is that G is the Green's function of the differential operator d b n . D = al - 6 V ' ( 4 - 2 ) The solution to the equation 4.1 is well known (for example, see reference [57]): e-3(r-r>)2/2b(s-s>) G ( r ' r ' S ' S ' ) = ( 2 ^ - ^ / 3 ) ^ ( 4 3 ) where d is the dimension of the space in which the polymer is embedded. 1 Recall that this is a polymer whose persistence length is much smaller than its contour length. 2 This refers to the position of the polymer a distance s from a specified end as measured along the contour 4. Single Chain Tethered to an Interface 52 In fact there are many different solutions, as any solution to the homoge-neous equation can be added to G to yield a new solution. In other words, if DH = 0, then G' = G + H is also a solution to the diffusion equation (4.1). Boundary conditions are required in order to specify a unique solution. To describe a polymer in the presence of an impenetrable plane, let us use cartesian coordinates in three-dimensions, with the barrier given by the xy-plane. Then we can represent r by (x,y,z) with the z-coordinate repre-senting the "height" above the barrier, and represent r' by the corresponding primed coordinates. Then the boundary condition is given by G(z = 0) = 0 [17]. As suggested in [17] and [58], this problem can be solved by the method of images. Just as in electromagnetism, we place a virtual source, or "image charge", of the opposite sign on the opposite side of the plane to enforce the required boundary condition. This yields: -3(r-r')2/2b(s-s') _ „-3(r-f') 2/2b(s-s') G < " > • • s '» = ( a r t ( . - " y ) / 8 ) V » <") where r' = (x',y', —z'). Height distribution We are primarily interested in the thickness and com-pressibility of the polymer loop, so integrating over x and y yields the prob-ability density for height: -3{z-z')2/2b{s-s') _ -3 (z+z ' ) 2 /26(s-s ' ) G ( " ' S'» " C M O - W * ( 4 ' 5 ) This result is very useful. For example, it allows us to calculate the normalized height density from the joint probability density, as presented in [7]. Let L be the contour length of the polymer, and 5b the height of the polymer ends. Then for a polymer loop with ends tethered infinitesimally close to the plane, , , ,. S^G(z,8b,s,0)G(5b,z,L-S,0)d8 viz) = hm ^ (4.0) S h ^ ° /o°°/o s, 0)G(5b, C,L-s, 0)dsdC p(£)d£ represents the ensemble average of the fraction of the polymer which lies between the planes z = Q and z = C + dC Carrying out the integrals in equation 4.6 yields: p(z) = l2{z/Lb)e-&z2lLb (4.7) which is simply a Rayleigh distribution. 4. Single Chain Tethered to an Interface 53 The moments of this distribution about zero are given by [7]: roo {zn) = / znp(z)dz (4.8) Jo ( L 6 / 6 ) n / 2 n ! ! ( 7 r / 2 "+ 1 ) 1 /2 n G dd (L6/6) n/ 2(n/2)! n even { ' The distribution can also be characterized by the descriptive statistics of the mean (fi), standard deviation (cr), skewness (c/i), and kurtosis (52) ([59, 60]). The mean is given by (z) as usual. The other three measures shall be defined in terms of the moments about the mean: W j = ( ( z - M ) j ) (4.10) Then we have: (4.11) (4.12) (4.13) o measures the spread or variance of the distribution about the mean. g\ measures the asymmetry of the distribution about the mean. It is simply the third moment about the mean, scaled by the standard deviation so that it is a dimensionless quantity, g^ is a measure of the flatness or peakedness of the distribution. Like g\, it is dimensionless due to scaling by the stan-dard deviation. 3 has been subtracted so that the kurtosis of the normal distribution is zero (this is fairly standard, but not always the case). Thus, if g2 is positive, then the distribution is more peaked than the normal distri-bution, and has higher tails (leptokurtic); if #2 is negative, the distribution has a flatter top than the normal distribution (platykurtic). For the normalized height distribution given by equation 4.7, these statis-0 = Vu2 91 = 1 ""2VW2 U4 92 = " 2 4. Single Chain Tethered to an Interface 54 tics evaluate to: 9i 92 = (2TT-6)V^ (4 - 7 r ) 3 / 2 0.6311 -6TT 2 + 24TT - 16 (4 - TT)2 0.2451 (4.14) (4.15) (4.16) (4.17) (4.18) (4.19) 4.1.2 The computer model The computer model that was used to study a single polymer loop was a modification of the model described in Chapter 3. I examined two particular parameter sets: the phantom chain with an intervertex bond stiffness of 100,. referred to as PPC100, and the corresponding self-avoiding chain, SA100. To prevent the polymer from passing through the plane, an excluded volume interaction was introduced between the spheres and the plane. The potential was an infinite square barrier when represented as a function of the separation distance; it was the same as the intervertex potential introduced in the self-avoiding chain. Thus, for both P P C and SA chains, the vertices acted as hard spheres of radius 0.42045 with respect to the barrier. The vertices at the two ends of the polymer were constrained to remain a distance 5b from the plane. I chose 5b to equal the excluded volume radius, 0.42045, so that the polymer was able to lie flat. To ensure that the end points did not change their height, the z-component of each trial move for these vertices was set to zero. In regard to the coordinates, the impenetrable plane was chosen to lie parallel to the xy-plane, at z = —0.42024. The end points remained in the xy-plane (z = 0) and all heights were measured with respect to these two vertices. As before, the first 10 configuration that were sampled were ignored, in order to allow the system to relax. Ensemble averages were generally calculated by sampling 200 post relaxation configurations for each parameter setting. One exception was in the determination of the height distribution for short polymers (N < 15), where 1000 configurations were used. 4. Single Chain Tethered to an Interface 55 4.2 Shape Properties of Loops as a Function of N I investigated the shape properties of a single polymer loop as a function of its contour length with the computer model. In particular, power laws were determined for a variety of measurements as a function of the number of segments, N. These results were compared with the results of Chapter 3, in order to determine the effect of tethering on the shape of the polymer. In addition, results were compared with the Diffusion Model presented in the previous section, by replacing the length parameter with the number of segments, that is L i—> Nb, where the segment length b is the unit of length. (r 2 e) and (r 2) were examined in order to determine the mean global length scale of the polymer. The mean end-to-end distance, (ree), was also measured, as it should represent the stress-free end-to-end distance of the polymer loop. Particular attention was paid to the thickness and height density of the polymer loop. For example, the mean top height ((top height)) was exam-ined, which was defined to be the height of the highest polymer segment, averaged over all samples. The variance of the top heights, a 2 o p h e i g h t was also measured, as together, these two measurements should give an approx-imation of the elastic compressibility of the polymer loop. In addition, the height distribution of the polymer loop was computed for different values of N, and the results compared to the analytical model. 4.2.1 Results and Analysis of Computer Simulations Simulations were performed for chains varying in length from 3 to 90 seg-ments. I was primarily interested in scaling laws that were indicative of large N behaviour; thus, power law fits were generally performed for N > N\, where TVi was chosen so that the log-log graph of the data looked linear within this regime. Due to the computational constraints on the size of N, the calculated scaling laws were more suggestive of universal behaviour than definitive. However, whenever the behaviour differed significantly from the expected theoretical values, additional simulations were performed for larger TV* to confirm or reject the power law. In figure 4.1 the characteristic length scales (r 2 e) and (r 2) are presented for both phantom and self-avoiding chains (PPC100 and SA100). A sum-mary of power laws fits is presented in Table 4.1. (r 2 e) scaled as expected, with the powers very close to the theoretical values for a free chain, as pre-sented in Chapter 3, table 3.2. In fact, they were equal within the stated 4. Single Chain Tethered to an Interface 56 Fig. 4.1: (a) (r2ee) and (b) (r2g) vs. number of segments (N) for P P C 1 0 0 and SA100 polymer loops. kc = 0. The data is presented in a log-log plot. A fit to the equation y = axb was performed for N > 15, and the results superimposed. errors. However, the prefactors were slightly smaller than both the theoret-ical and computed values given in Chapter 3. In other words, the interface appeared to bring the ends closer together on average, but always by the same factor. On the other hand, for the radius of gyration, (r 2), both the exponents and the prefactors were quite close to the values expected for a free chain in solution. It is also interesting to note that the power law appeared to be a reasonable fit to the data, all the way down to JV = 4. The linear measures of mean end-to-end distance and mean top height are presented in Figure 4.2, with the results of the power law fit in Table 4.1. The fit was performed for ./V between 15 and 90. As expected, the scaling exponent for (r e e) was half the size of the (r 2 e) exponent, within error. The mean top height grew more quickly with N, but with a smaller prefactor. Fluctuations in the top height are interesting as they can provide information about the thickness compliance of the polymer loop. In particular, if the distribution is Gaussian, then the top height can be roughly represented by an effective Hookean spring. This allows the effective spring constant, or 4. Single Chain Tethered to an Interface 57 Mean value PPC100 SA100 Figure of observable a b a b r2 ' ee 0.68(8) 0.99(3) 1.0(1) 1.17(3) 4.1(a) r2 9 0.17(1) 0.95(2) 0.172(8) 1.20(1) 4.1(b) ree 0.74(5) 0.49(2) 0.92(5) 0.59(2) 4.2(a) top height 0.47(2) 0.58(1) 0.39(1) 0.72(1) 4.2(b) 2 '"'top height 0.024(5) 0.98(6) 0.015(3) 1.26(5) 4.3(b) 19(2) -0.40(4) 26(2) -0.55(2) 4.4 Tab. 4.1: Power law fits for P P C 1 0 0 and SA100 single polymer loops. The N dependence of various quantities was investigated by fitting the data to y = aNb for N between 15 and 90. Number of segments Number of segments Fig. 4.2: (a) ( r e e ) and (b) (top height) vs. number of segments (N) for P P C 1 0 0 and SA100 polymer loops. kc = 0. The data is presented in a log-log plot. A fit to the equation y superimposed. ax was performed for N > 15, and the results 4. Single Chain Tethered to an Interface 58 the elastic modulus, to be determined from fluctuations in the top height. I defined a top height elastic modulus, Mj_, as: M X = kDT{t)/({t2) - (t)2) (4.20) where t is the top height, and T is the temperature. In our system of units, ksT is unity, so M± is simply the reciprocal of the variance in the top height scaled by the mean top height. In the simulation, the fluctua-tions in the top height were approximately Gaussian, as presented in Figure 4.3. The variance in the top height was examined as a function of N. Er-rors were estimated by the standard error of the variance, which is simply \/2 x variance/VSample size [59]. (Note: this formula was also used to estimate the error in the denominator of M±.) The data was fit to a power law which appeared to be valid within the 15 < N < 90 regime. The results are in Table 4.1. A plot of Mj_ is presented in Figure 4.4. Results of a power law fit also appear in Table 4.1. Next, the normalized height densities were examined. A few character-istic normalized histograms (frequency polygons) for PPC100 are presented in Figure 4.5, and for SA100 in Figure 4.6. As expected, the height ap-proximately follows a Rayleigh distribution, when the number of segments was sufficiently large; eight segments appeared to be enough. The distribu-tions were characterized by the four descriptive statistics mentioned earlier. Graphs of the mean height, standard deviation, skewness and kurtosis for PPC100 and SA100 are presented in Figure 4.7. Predictions from the Dif-fusion Model are superimposed. From these figures the distributions appear to approach stable values for N > 8. Power law fits for the mean and standard deviation for N > 8 are summarized in Table 4.2. In addition the mean values for g\ and gi within this range are also given. The mean height of the PPC100 and SA100 polymers grew approxi-mately 10% more quickly than (ree). This was likely due to the steric re-pulsion effectively pushing the polymer away from the plane barrier. The standard deviations grew in a similar manner, but not quite as quickly with N. One surprise was that the phantom polymer did not appear to reach ideal behaviour by N = 90 for the mean height. Instead, the power law was larger than the predicted value of 0.5. This was expected for short chains, as the height must go to zero when the number of segments is 1, but it was surprising that unlike most of the other measurements, such as the standard deviation of the height, 90 segments was still not enough. According to the figure, the intersection of the PPC100 and D M power laws occurs for N 4. Single Chain Tethered to an Interface 59 Top height Number of segments Fig. 4.3: Variance of top height ( < 7 t o p height ) for P P C 1 0 0 and SA100 polymer loops. kc = 0. (a) Histogram of top height measurements for a P P C 1 0 0 loop wi th N — 90. A normal curve that was determined by the mean and variance of the data is superimposed, (b) ot0p height vs. number of segments (N), presented in a log-log plot. A fit to the equation y = axb was performed for > 15, and the results superimposed. Mean value PPC100 SA100 Figure of observable a b a b Mean 0.270(4) 0.552(4) .248(5) 0.664(5) 4.7(a) Std Dev 0.173(3) 0.520(5) 0.143(4) 0.665(8) 4.7(b) Skewness (gi) 0.63(2) - 0.54(1) - 4.7(c) Kurtosis (52) 0.16(6) - -0.18(3) - 4.7(d) Fitting parm., Pi 3 0.58(3) 1.10(1) 0.48(2) 1.34(1) 4.9 Tab. 4.2: Characterization of the height distr ibution of P P C 1 0 0 and SA100 single polymer loops. The N dependence of the first three descriptive statistics was investigated by fitting the data to y = aNb for N between 8 and 90. For the skewness and kurtosis, the table simply presents the mean value, calculated from this same range of N. 4. Single Chain Tethered to an Interface 60 1 -j 1 1 1—i—i i i i | 1 1 1—i—i i i i | 1 1 10 100 Number of segmen ts Fig. 4.4: Top height modulus ( M i ) vs. .number of segments (N) for PPC100 and SA100 polymer loops. kc = 0. A fit to the equation y = axb was performed for N > 10, and the results superimposed. 4. Single Chain Tethered to an Interface 61 Fig. 4.5: Height distribution for PPC100 polymer loops with kc = 0. (a) His-togram of the normalized probability density for N = 90 with a fit to the Rayleigh distribution superimposed, (b) Normalized histograms for the height distribution of chains of various lengths (JV 6 {4,15,45, 90}). Note the difference in vertical scale between graphs (a) and (b); the histogram for N = 90 is the same in both graphs. 4. Single Chain Tethered to an Interface 62 6 8 10 12 14 16 Height 2 4 6 8 10 12 14 16 Height Fig. 4.6: Height distribution for SA100 polymer loops with kc = 0. (a) Histogram of the normalized probability density for N = 90 with a fit to the Rayleigh distribution superimposed, (b) Normalized histograms for the height dis-tribution of chains of various lengths (N e {4,15,45,90}). Note the dif-ference in vertical scale between graphs (a) and (b); the histogram for N = 90 is the same in both graphs. 4. Single Chain Tethered to an Interface 63 Fig. 4.7: Descriptive statistics for the height distribution of polymer loops, as a function of segment number N. Polymer were of type PPC100 and SA100, with kc = 0. (a) Mean, (b) standard deviation, (c) skewness (gi), (d) kurtosis (32)- In graphs (a) and (b), a fit to the equation y = axb was per-formed for N > 8, and the results superimposed. In addition, predictions of the Diffusion Model were superimposed as dashed lines. 4. Single Chain Tethered to an Interface 64 0.1 -1 1—i—i—i i 1111 1—i—i—i i 1111 i—i—i—i i 1111 1 1 10 100 1000 Number of segments Fig. 4.8: Mean height vs. number of segments for a PPC100 polymer loop. Pre-dictions of the Diffusion Model and the power law determined from N between 8 and 90 are superimposed. in the low hundreds. It was expected that a transition to ideal behaviour would occur at approximately this number of segments. To investigate this possibility, additional simulations were performed with the PPC100 model for N ranging from 200 to 800. Similar simulations were not performed for the self-avoiding polymer due to limited computer time. These results are presented in Figure 4.8 together with the predic-tion of the Diffusion Model and the power law determined from the PPC100 model for N between 8 and 90. As expected, the data appeared to roughly follow the D M at large values of N, with the transition occuring at approx-imately Af = 200. However, in the range of N tested it was quite difficult to distinguish between the two different power laws, so the results are not definitive. Simulations of much longer loops would be required for a stronger confirmation of this behaviour. The mean value of the skewness (gi) was close to the theoretical value 4. Single Chain Tethered to an Interface 65 for the phantom chain within the 8 < TV < 90 range. On the other hand, the self-avoiding chain was slightly more symmetric. As the number of segments decreased below 10, both the PPC100 and SA100 distributions became increasingly symmetric. The kurtosis measurements indicated that the phantom chain was slightly more flat-topped than predicted by the Dif-fusion Model, with the self-avoiding chain even more flat-topped. These deviations from the analytical model became more pronounced as the num-ber of segments decreased below 10. In addition to these statistical measures, the numerical height distribu-tions were fit directly to the Rayleigh distribution (Equation 4.7) in the fitting parameter P\ = Lb for A^ > 8. Recall that L represents the con-tour length of the D M loop, and b the mean step size. Thus, P\ was called the "effective Nb2". The results of the fit are presented in Figure 4.9. P i appeared to follow a power law in N, but grew faster than the 1 : 1 cor-respondence predicted by the Diffusion Model. However, the results were consistent with the mean height data, as the exponents for the effective Nb2 were approximately twice as large as the exponents for the mean height. 4. Single Chain Tethered to an Interface 66 1 10 100 Number of segments Fig. 4.9: Effective Nb2 (P,) vs. number of segments (AT) for PPC100 and SA100 polymer loops. Values were determined by fitting the normalized probabil-ity densities to the Rayleigh distribution v(z) = f2(z/Pi)exp(—Qz 2/P\). The data is presented in a log-log plot. A fit to the equation y = axb was performed for N > 8, and the results superimposed. In addition, the predictions of the Diffusion Model are also superimposed as a dotted line. 4. Single Chain Tethered to an Interface 67 4.3 Shape Properties of Loops as a Function of the Bending Stiffness In this section, I studied the shape properties of polymer loops and their dependence on bending stiffness. This included a phenomenological exam-ination of the appropriateness of representing a worm-like chain by i) the diffusion model, and ii) an equivalent freely jointed chain. Due to the nature of the data, this discussion is more qualitative than the previous section. As described in the previous chapter, a W L C can be represented by an equivalent F J C chain consisting of N' statistical segments of length 6', pro-vided the contour length L is sufficiently large with respect to the persistence length lp. b' is known as the Kuhn length, and for a free W L C , b1 = 2lp. Furthermore, in section 3.3, a relationship was found between the compu-tational parameter kc, and the equivalent bending stiffness coefficient keg, where fceff = lp (within our system of units), just like the bending stiffness coefficient of the ideal W L C . 4 These correspondences allow us to transform the results of the previous sections into predictions for the present section. In particular, an observable given as a function of N can be re-written as a function of kc. Let b' = ab where a is a dimensionless quantity, and b is the segment length of a F J C or D M chain with no explicit bending stiffness. Now the Kuhn equivalent chain has the same length as the original chain, so L is unchanged, and N' = N/a. Then given a function of N and b for a F J C (or D M chain) with no bending stiffness, we obtain a result for the Kuhn statistical chain by making the substitutions: b i — > ab N i—> N/a (4.21) For example, consider the height distribution of the Diffusion Model. Letting a = 2kes yields l> = \ J J V N £ H (4.22) a _ J ^ - ^ (4 2 3 ) from equation 4.14. Note that b is not explicitly present, as in this system of units, b is effectively 1. As for the skewness and kurtosis, they are unaffected by the scaling. 4 Specifically, it was found that to a reasonable approximation, fceff = 0.85(9)+ 1.09(4) x kc for the PPC100 phantom chain, and fceff = 1.0(1) + 1.12(4) x kc for the SA100 self-avoiding chain. 4. Single Chain Tethered to an Interface 68 As another example, let us consider the power law relationship: (M) = aNpbq (4.24) where M has units of length to the power of q. Then in terms of kes, in our natural system of units. Now, let us consider an added complication: the existence of excluded volume. Naively, one could approximate a W L C with excluded volume by an excluded volume FJC using the substitutions given above. However, this fails to take into account the implicit bending stiffness which results from the steric repulsion between segments in an excluded volume FJC; the effective (statistical) segment length of this chain is larger than simply b. To correct for this, a in equation 4.21 should be divided by the Kuhn length of the excluded volume FJC. For concreteness, let lo be the persistence length of a self-avoiding F J C with no explicit bending stiffness (i.e. kc = 0). Suppose that equation 4.24 is a scaling law for this chain. Then equation 4.25 still represents an approximation scaling law for a W L C , provided we make the substitution ketf i—> keg/(2lo). I will refer to this scaling as SA SS, or self-avoiding "smart scaling". The IQ which was used in calculations came directly from the angular correlation exponential fit of Chapter 3, yielding J 0 = 0.82(3). However, it is expected that when kc is sufficiently large, but still much smaller than the contour length, the excluded volume W L C can be accu-rately represented as phantom FJC . As mentioned in Chapter 3, at a given contour length, when the persistence length is much larger than the size of the excluded volume, the polymer behaves effectively like a phantom chain. It is important to note that all of these scaling relations must fail when kc/L is sufficiently large, as the W L C will approach the behaviour of a rigid rod. Furthermore, the relationships between fceff, lp and the Kuhn length were all established for a free W L C in solution, and it is not clear that they will hold under the tethered boundary conditions. It is interesting to examine when this approximation is valid, and when it breaks down. It is important as well, since the F J C is much easier to simulate than the W L C . Simulations were performed for A^ = 90, and kc ranging from 0 to 10, for both the phantom and self-avoiding chains. Various quantities were mea-sured, and plotted against the effective bending stiffness coefficient fceff. (M) = (2keft)"-PaNP (4.25) 4.3.1 Results and Analysis of Computer Simulations 4. Single Chain Tethered to an Interface 69 350-, 0 2 4 6 8 10 12 0 2 4 6 8 10 12 Effective bending stiffness coefficient Effective bending stiffness coefficient Fig. 4.10: (a) (rf;e) and (b) (r2g) as a function of the effective bending stiffness coef-ficient, fcoff, for P P C 1 0 0 and SA100 polymer loops. N = 90. Predictions of the P P C , S A , and S A SS scaling models have been superimposed. Figure 4.10 shows the global length scales (r 2 e) and (r2) as a function of the bending stiffness. The power laws determined from the simulations in section 4.2 were converted into functions of fceff using the scaling arguments given above, and these predictions superimposed on the data. As expected, the equivalent phantom F J C provided the best fit, closely following (r 2 e) for all values of keff. On the other hand, the naive SA scaling predictions deviated almost immediately from the data. Rescaling by the persistence length of the SA chain improved the fit, as seen from the SA SS predictions. For the (r2) data, P P C scaling yielded the best fit, as before. However, the data appeared to deviate from the predicted behaviour at around keg = 4, corresponding to a Kuhn lenght of 8, and hence approximately 11 statistical segments. The linear measures of the mean end-to-end distance ({ree)) a n d the mean top height are presented in Figure 4.11. P P C scaling appeared to yield the best fit for (ree) for both the PPC100 and SA100 data. For the mean top height, this was also true over the domain fceff < 4. However, above this stiffness, the polymer height grew more slowly than predicted, eventually meeting up with SA SS predictions. For kes > 8 (corresponding to 4. Single Chain Tethered to an Interface 70 <H—•—r~'—i—1—i—•—i—•—i—1—i—1 °H— 1 —i— 1 —i—•—i— 1 —i— 1 —i— 1 —i— 1 0 2 4 6 8 10 12 0 2 4 6 8 10 12 Effective bending stiffness coefficient Effective bending stiffness coefficient Fig. 4.11: (ree) and (b) {top height) as a function of the effective bending stiffness coefficient, fccff, for P P C 1 0 0 and SA100 polymer loops. N = 90. Predic-tions of the P P C and S A SS scaling models have been superimposed. 5 or 6 statistical segments), the mean top height began to decrease, perhaps signifying the start of a transition to rigid rod behaviour. However, this was not definitive due to fluctuations in the data. It is also interesting to note, that for keR > 3, it was difficult to distinguish between the phantom and self-avoiding data. This was consistent with our predictions for SA polymer behaviour when the persistence length was much larger than the excluded volume sphere radius (i.e. lp ^> Rs). Next, the variance of the top height was examined (Figure 4.12), to gain some insight into how bending stiffness affects the elastic compressibility. Similar to the mean top height, the SA SS predictions appeared to match the data most closely, for heg > 3. A plot of the top height modulus, Mj_ appears in Figure 4.13. Characteristic plots of the height densities for different values of kc are presented in Figure 4.14. Consistent with the changing N data, for small values of kc, the SA distributions were wider and flatter than the P P C distributions. As kc became larger, the P P C distributions also began to look wider and flatter. In addition, there was a smaller transition from 4. Single Chain Tethered to an Interface 71 Fig. 4.12: Variance of top height as a function of the effective bending stiffness coef-ficient, fcoff, for PPC100 and SA100 polymer loops. N = 90. Predictions of the PPC and SA SS scaling models have been superimposed. 4. Single Chain Tethered to an Interface 72 PPC100 SA100 P P C scaling SA SS T 2 4 6 T 8 1— 10 12 0 Effective bending stiffness coefficient Fig. 4.13: Top height modulus ( M x ) as a function of the effective bending stiffness coefficient, fcoff, for P P C 1 0 0 and SA100 polymer loops. N = 90. 4. Single Chain Tethered to an Interface 73 (a) PPC (b) SA Height Height Fig. 4.14: Height distributions for (a) P P C 1 0 0 and (b) SA100 polymer loops wi th various settings of kc. N = 90. Normalized frequency polygons are shown for kc £ {0 ,1 ,2 ,6} . kc = 0 to kc = 1 for SA than for P P C . This was likely due to the inherent bending stiffness in the self-avoiding chain resulting from steric repulsion. Changes in the height distributions were quantified by examining the de-scriptive statistics as presented in Figure 4.15. Consistent with the changing N results, the mean height never really agreed with the Diffusion Model, even for kc = 0. Consistent with the previous graphs, P P C scaling matches the data most closely for keg between 1 and 4. This was true for both the PPC100 and the SA100 data. Above fceg- s=s 4 for the phantom chain, and fceff « 5 for the self-avoiding chain, the data deviates from the P P C scaling curve, growing more slowly than predicted. Above these values, the SA SS scaling prediction bends closer to the data points. The standard deviation of the height followed the same basic story. However, the primary difference was that the Diffusion Model predictions were actually reasonably good for a up to A;eff = 4. Both the skewness and the kurtosis remained close to the kc = 0 values for small ken, then fell below these values for larger ke^, consistent with the small chain behaviour. To further investigate the deviation from the equivalent phantom F J C behaviour, another set of simulations was performed for a phantom W L C 4. Single Chain Tethered to an Interface 74 —\— 10 0 2 4 6 Effective bending stiffness coefficient (b) (c) • A -* • - —•—PPC100 —A— SA100 • . i i i i . DM i • i ' i 2 4 6 8 10 12 Effective bending stiffness coefficient -PPC100 h-A— SA100 ---- DM — PPC scaling SA SS 0 2 4 6 8 10 12 Effective bending stiffness coefficient 0 2 4 6 8 10 12 Effective bending stiffness coefficient Fig. 4.15: Descriptive statistics for the height distr ibution of P P C 1 0 0 and SA100 polymer loops as a function of the effective bending stiffness coefficient, fceff. N = 90. Predictions of the Diffusion Mode l , and the P P C and S A SS scaling models, have been superimposed. 4. Single Chain Tethered to an Interface 75 of length N = 45. The results of this study are presented in Figure 4.16. For comparison, the results of the N = 90 W L C were scaled by a constant factor so that they would agree at kc = 0. These results are superimposed, as are the predictions of P P C scaling and D M scaling. Overall, the predictions of P P C scaling and the D M model were not as good for the N = 45 chain as for the N = 90 chain. This was not surprising as these theories should primarily describe large N behaviour. As before, the predictions of the D M model were better for the standard deviation than for the mean. It was expected that the N = 45 chain would deviate from the P P C scaling behaviour at a fceff value approximately half of the corresponding fceff value for the longer chain. This appeared to be roughly true, although there was not enough data to determine this with confidence. In addition, it should be noted that while the short chain data followed the basic curvature of the P P C scaling results for fceff < 3, the fact that the two plots did not agree even for kc = 0 made it difficult to compare them directly. As for the skewness and kurtosis data, the short chains appeared to decrease from the D M values slightly more quickly than the long chains, but due to noise in the data, this too was difficult to evaluate. So to summarize, the equivalent F J C appeared to be a reasonable rep-resentation of a W L C polymer loop when there were approximately 11 or more statistical segments. On the other hand, while the Diffusion Model was able to predict the basic shape of the height density quite well, it was not particularly good for the mean height when N was less than 90. It did not appear that we were quite in the large TV regime with respect to this measure. 4. Single Chain Tethered to an Interface 76 Fig. 4.16: Descriptive statistics for the height distribution of a P P C 1 0 0 polymer loop as a function of the effective bending stiffness coefficient, fceff. N = 45. The results for N = 90 have been superimposed, but scaled to agree wi th the N = 45 chain for kc = 0. Predictions of the Diffusion M o d e l and the P P C scaling model have been superimposed. 5. Single Loop under Parallel Confinement 77 5. S I N G L E L O O P U N D E R P A R A L L E L C O N F I N E M E N T The focus of this chapter is the elastic compressibility and shape character-istics of a single polymer loop under parallel confinement. Various models of confined polymers are presented in the first section. The results of the computer simulation and comparisons with analytical predictions are in the second section. 5.1 Models and Theory The specific system of interest is a single polymer confined between two in-finite planes, with both ends tethered to a single plane. As before, the ends are free to slide along the plane. In the present section, the polymer models introduced in the previous chapters are applied to this system. Analytical predictions from the Diffusion Model are discussed first, followed by stan-dard results for untethered polymers. A discussion of how the system was modelled computationally concludes the section. 5.1.1 The Diffusion Model Recall the Diffusion Model introduced in Chapter 4. The probability dis-tribution (or Green's function) for this model must satisfy the differential equation 4.1 subject to the appropriate boundary conditions. In the case of a polymer confined between two planes, we must not only enforce the condition G(z = 0) = 0, but also the condition G(z = h) = 0, where h is the gap distance. As before, a solution may be found by using the method of images. However, in this case the resulting solution is an infinite sum, due to the "infinite reflections" of the image charges in the two facing planes. The solution is given by: CO G{z,z',s,0;h) = ^/2,/2-KLb ^ {exp[-3(z - z' - 2nh)2/2Lb] n=—oo - exp[-3(z + z' - 2nh)2/2Lb}} (5.1) 5. Single Loop under Parallel Confinement 78 where L is the contour length of the polymer and b is the effective segment length as before [17, 58]. This series converges rapidly when the separation distance h is large (h 2> \/Lb). An equivalent series solution that converges quickly for small separations (h <C y/Zb) is: CO G(z, z', s, 0; h) = (2/h) ^ s'm(mrz/h) sm(nnz'/h) 71=1 x exp{-Lbn2n2/6h2} (5.2) This solution may be found through the use of an eigenfunction expansion [17, 61]. These results are discussed in Heinrich et al. [7] and used to calculated the free energy of compression for a flexible polymer loop. I will summarize the calculation here. Consider a polymer loop of length L. Let f2(L, h) be the number of configuration of the loop confined within a gap distance of h, relative to the number of configurations of the unconfined loop. That is, fi(L, ti) = lim G{5b, 5b, L, 0; h)/G{5b, Sb, L, 0; oo) (5.3) where G(8b, 5b, L, 0; oo) is the probability density for the unconfined loop, as calculated in Chapter 4 (equation 4.5). Then the change in entropy due to confinement is simply: A S = ks log(fi(L, ti)) (5.4) As there is no internal energy in the Diffusion Model, the free energy depends only on the entropy. Thus, the free energy of compression for a single loop is given by: e = -kBTlog{n(L,ti)) (5.5) The expressions for G in equations 5.1 and 5.2 can be used to approx-imate £l(L,h). Replacing each infinite series by the lead term, then ex-panding to quadratic order in 5b yields a good approximation within the appropriate regime. Calculating the free energy of compression from this approximation gives: eS{h) « kBT[ir2Lb/6h2 + 3 ]n(h/y/Tb) - ln(27r5/27)/2] (5.6) in the small h regime, and eL(h) « - jfeBT{ln[l - (24h2/Lb - 2) exp(-6/i2/Lb)}} (5.7) 5. Single Loop under Parallel Confinement 79 <H 1 , — , 1 , 1 OH , 1 , r ^ - < 1 1 1 0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0 h / ( L b ) " 2 h / (Lb ) 1 ' 2 Fig. 5.1: (a) Free energy and (b) force vs. gap distance for the Diffusion Model. Curves for the large and small gap distance regimes have been superim-posed. Values have been scaled appropriately by \/~Lb and ksT to yield dimensionless quantities. in the large h regime. These two expressions for the energy are plotted in Figure 5.1(a). While the curves do not intersect, in the region h/VTb G [0.70,0.76] they differ by less than 0.001 kBT. \es{h) - eL(h)\ is a minimum at gap distance of hc « OJSVLb-To hold the two planes at a given gap distance h, an inward compressive force, F, is required on each plane. F can be calculated from the free energy of compression by taking the partial derivative with respect to the gap distance: F = —de/dh. The expressions given above for e(h) yield: Fs(h) » kBT[TT2Lb/3h3 - 3/h] (5.8) for h \fh~b, and Fr (h) ~ { h / L b ) eM-Gh2/Lb)(4h2/Lb - 1) L [ ' ~ 24{h?/Lb) exp(-6/ l 2 /L6) - 2exp(-6/ l 2 /L6) - 1 [ ' for h >^ VL6. A plot of these force curves is presented in Figure 5.1(b). As the curves intersect when the gap distance is hc, this is an appropriate 5. Single Loop under Parallel Confinement 80 switch-over point for a piecewise construction of F(h) valid for all h > 0. That is, I shall approximate the compressive force F(h) by Fs(h) when h < hc and by Fr{h) when h > hc. The energy function e(h) for h > 0 may be constructed in the same way. When the gap distance is small, the compressive energy and force should be the same as for a flexible untethered polymer. This will be shown ex-plicitly in the next section, which investigates untethered polymers under parallel confinement. Many properties of confined untethered polymers have been calculated ana-lytically that have yet to be established for polymer loops. In particular, the effects of bending stiffness and excluded volume, which are not addressed by the Diffusion Model, have been studied in certain regimes. However, let us first examine the differences between tethered and untethered polymers under confinement. Using Green's functions, the partition function of a flexible polymer confined in a box has been calculated [16, 18]. This allows the free energy of compression and the corresponding compressive force to be determined in the standard way, as discussed for the Diffusion Model. In the limiting case when one dimension of the box is much smaller than the characteristic size of the polymer (\/Lb), and the other two dimensions are much larger than the polymer, the force on the confining planes is given which is the same as the leading term in equation 5.8 for the tethered poly-mer. When h is sufficiently large, the relationship between compressive force and gap distance reduces to that of an ideal gas. However, in the tethered polymer Diffusion Model, the force falls off much more quickly due to the Gaussian factor which dominates for large h, as evident from equation 5.9. The behaviour of a confined untethered polymer changes significantly when the bending stiffness is increased but the contour length is held fixed. When the polymer is stiff, the relationship between force and gap distance can be calculated by considering the exclusion of bending modes due to lateral confinement [62]. This results in the scaling law: 5.1.2 Untethered polymer models by: F{h) = kBTir2Lb/3h3 (5.10) F{h) ~ knT/h5/3 (5.11) 5. Single Loop under Parallel Confinement 81 It is not clear how well this description holds for tethered polymers, even under strong confinement. However, the cross-over from the flexible to rigid regime should also be manifest in the relationship between F and h. Another description which may be used to study the free energy and shape properties of confined polymers is the "blob model" of de Gennes [10]. While this simple scaling model applies only to flexible chains, it is very useful because it can take into account the effects of excluded volume. Given a linear chain with N segments, the basic idea is to divide it into a string of N/rj blobs containing rj segments each. Then within each blob, the effects of the boundaries are to be treated as a weak. This will be illustrated with an example. Consider a phantom chain consisting of N segments of size b, confined between two parallel planes with a gap distance of h. Let us break the polymer into blobs of diameter h. Assuming that within each blob the effects of confinement are weak, h ~ \frjb, where rj is the number of segments within each blob. (For this assumption to hold, h 3> b.) Then (r2e) for the total polymer is given approximately by (N/rj)h2 ~ Nb2. As expected, (r2e) for a flexible phantom chain is not affected by parallel confinement, as long as the gap distance is much larger than the segments size [10]. This is essentially because the power scaling law for (r2e) is the same in two dimensions as in three. Things are more interesting for the self-avoiding chain. Using the ar-guments given above, but the Flory expression 3.14 for the characterstic global size, we obtain the relationships: h ~ 7 ? 3 / 5 & , and (r2e) = (N/n)3/2h2. Combining these expressions yields the scaling law: ^ 3 / 4 5/4 (d) ~ ^ f - (5-12) The blob model can also be used to estimate the confinement energy of a flexible self-avoiding chain [10]. As argued by de Gennes, in the strong confinement regime, this gives the scaling law: e{h) ~ kGTN(b/hflz (5.13) Hence, the force is given by: F(h) ~ kBTNb5/3/h8/3 (5.14) 5.1.3 The computer model To simulate a confined polymer loop, a second infinitely hard plane was added to the computational model. As before, the vertices interacted with 5. Single Loop under Parallel Confinement 82 the plane through an infinite square well potential (i.e. the chain vertices could not come closer than one sphere radius to the plane). The plane was allowed to move "up and down" along the z-axis, adding one additional degree of freedom to the model. It was treated as another body in the simulations, but with a unique reduced step parameter, D T B , which was set to 0.0005. This was ten times smaller than the DT parameter setting for the vertices. The moving plane was called the upper plane, and the gap distance was equated with the height of this plane. With respect to coordinates, it should be noted that the lower tethering plane was set to a height of z = —0.42045 as before. Thus, the minimum possible height of the upper plane due to sterics was z = 0.42045; this resulted in a small offset. Most of the simulations were performed under constant force (CF). That is, the upper plane was placed in a potential that increased linearly with height. The position of the upper plane was expected to fluctuate, as it was treated like another body in the simulation, ft was anticipated that the fluctuations would be quite small in the strongly confined (highly com-pressed) regime, and relatively large in the weakly confined regime. As I was primarily interested in the strongly confined regime, most of the simu-lations were performed under constant force. In order to examine the force when the system was under weak confinement, a slightly different simulation approach was used. The upper plane was placed in a harmonic potential to reduce the fluctuations in the position (HF); this modelled an upper wall connected to a spring. The price that was paid was a corresponding fluctu-ation in the force. Thus, by tuning the spring constant of the upper wall, fluctuations between force and position could be traded off of each other. PPC100 phantom chains and SA100 self-avoiding chains were used in all simulations. 5.2 Results and Analysis of Computer Simulations First, I will examine the shape distortions of a single polymer loop under confinement. Next, the dependence of the compressibility on the number of segments and the bending stiffness will be investigated. 5.2.1 Shape distortions and compressibility of a long flexible loop In this subsection, PPC100 and SA100 polymer loops containing 90 segments were confined under both a constant force (CF) and a harmonic force (HF), and the shape properties examined. For the C F simulation, the gap distance was varied from 1 to 12 (bond length units). No explicit bending stiffness 5. Single Loop under Parallel Confinement 83 was used (A;c = 0). The investigation focused on the effect of compression on the thickness and global size dimensions of the loops. Various measures of the thickness of the polymer loop as a function of force are presented in Figure 5.2. The force is given in scaled units of kBT/b. Data is presented for both the constant and harmonic compressive force sim-ulations. At all levels of force, the average top height was approximately twice that of the average mean height. In the strong confinement regime, the mean height of the upper plane (i.e. the mean gap distance) accurately tracked the top height of the polymer, separated by roughly a standard de-viation in the top height. Both the CF and the HF simulations exhibited this behaviour. This regime corresponded to a gap distance less than ap-proximately the mean unconfined top height of the polymer, or alternatively to forces greater than approximately unity. Furthermore, within this range of values, there was good agreement between the PPC100 simulations and the predictions of the Diffusion Model. However, in the low force/weak confinement regime, the C F simulations followed neither the Diffusion Model nor the rough shape of the mean top height. This was likely due to the large fluctuations in the gap distance. As seen in Figure 5.2(a,b), the mean top height and the mean gap distance were separated by approximately one standard deviation in the gap distance. This resulted in a roughly power law dependence between the force and the mean gap distance even for low values of force. Thus, for increasing gap distance, the force did not fall off as quickly as predicted by the Diffusion Model. In this regime, the HF simulations provided better tracking of the top height at the expense of an error in the force. At the largest values of gap distance, the relative errors in the forces were quite large due to fluctuations. Overall though, the PPC100 HF simulation followed the top height and the D M predictions quite accurately. As I was interested in the compressive properties in the strongly confined regime, a power law fit to the constant force data was performed for F > 0.7. For the PPC100 chain, the scaling law was F = aNb with a = 261(28) and b = 3.01(5). As before, the number in brackets denotes the error in the least significant digit. The exponent matches the leading term of the D M results within error. The prefactor was approximately 10% smaller than the coefficient of the leading term in the D M model, 907r2/3 = 296. This small discrepancy may have resulted from the omission of higher order corrections. It is also possible that 90 segments were not enough to put the chain into the long, flexible regime described by the D M model. This would be consistent with the results from Chapter 4. For the SA100 chain, the parameters of the power law relationship were found to be a = 210(17) and b = 2.57(4). 5. Single Loop under Parallel Confinement 84 (a) P P C ( C F ) X 0.01 0.1 1 10 Compressive force (c) P P C (HF) I III 0.01 0.1 1 Compressive force (b) S. • Gap distance • Top height * Mean height DM gap distance 0.01 0.1 1 10 Compressive force (d) SA (HF) 0.1 1 10 Compressive force Fig. 5.2: Height measurements as a function of force for (a,c) PPC100 and (b,d) SA100 confined polymer loops compared with the Diffusion Model (shown by dashed curves DM). Vertical error bars represent the standard devia-tion of the heights. Errors in the force are represented by horizontal error bars in the gap distance data. Figures (a,b) show compression under con-stant force (CF); figures (c,d) show the upper plane in different harmonic potentials (HF). Gap distances from (a,b) are superimposed on (c,d) (re-spectively) as open squares to facilitate direct comparison. In (a,b), a fit to the equation y = axb was performed for the gap distance for F > 0.7. 5. Single Loop under Parallel Confinement 85 Fig. 5.3: Height distributions for (a) PPC100 and (b) SA100 confined polymer loops under different settings of constant, compressive force. N = 90. Normal-ized frequency polygons are shown for F € {0, 0.3,1, 3,10}. The exponent was slightly smaller than 8/3, the value predicted by the blob model for a flexible untethered polymer. The height distribution of a loop under force is presented in Figure 5.3. This provides a graphical representation of the thickness of a polymer loop under compression/confinement. Next, (r-ge) and ( r 2 ) were measured as a function of gap distance to assess the effects of confinement on the global shape dimensions. This is presented in Figure 5.4. As expected, ( r 2 e ) and ( r 2 ) were independent of the gap distance for the phantom PPC100 chain. The self-avoiding chain SA100 showed roughly power law dependence on the gap distance in the strong confinement regime as predicted by the blob model. A power law fit was performed for h < 9, to determine the approximate scaling law. For ( r 2 e ) — ahb, the parameters were found to be: a = 874(83) and -0.57(6). The exponent differed slightly from the value of —1/2 predicted by the blob model, but it was reasonable for a rough measurement. The fitting parameters for ( r 2 ) = chd were: c = 102(5) and d = -0.42(3). 5. Single Loop under Parallel Confinement 86 (a) • PPC100 A SA100 at -a (b) ****** *** ** • PPC100 A SA100 M e a n g a p d i s t a n c e M e a n g a p d i s t a n c e Fig. 5.4: (a) (r2ee) and (b) (r2g) vs. mean gap distance (h) for P P C 1 0 0 and SA100 polymer loops. N = 90 and kc = 0. The data is presented in a log-log plot. A fit to the equation y •• was performed on the SA100 data for (h) < 9, and the results superimposed. 5. Single Loop under Parallel Confinement 87 5.2.2 Dependence of compressibility on N and kc under strong confinement The relationship between the compressive force and the gap distance was investigated in the strongly confined regime for polymer loops of different lengths, and different values of the bending stiffness. The constant force ensemble (CF) was used for all simulations in this subsection. Dependence on N PPC100 and SA100 polymer loops were compressed under constant force for different values of N. The value of N was varied from 15 to 90. No explicit bending stiffness was used (kc = 0). A representative family of curves for the phantom chains is presented in Figure 5.5. The mean top height curves were approximately parallel for all values of force. This was also true for the mean gap distance in the strong confinement regime. However, for F < 0.7, the gap distances data is less regular, due to errors induced by the large fluctuations. The curves for the self-avoiding chain were similar. Power law fits were performed for F > 0.7, and the results summarized in Figure 5.6. The exponent was approximately independent of the number of segments within the range examined. Thus, N = 15 appeared to bring us into the large N regime with respect to the force-gap distance scaling law. The average value of the exponent was 3.00(2) for the phantom chain, and 2.64(2) for the self-avoiding chain. This was in close agreement with the predictions of the random flight (DM) and blob models, respectively. Thus, the effect of tethering on the force vs. gap distance exponent was small in the strongly confined regime, even for N as small as 15. The prefactors grew approximately linearly with /V as expected. While the values for the PPC100 chain were smaller than predicted by the Diffusion Model, this was consistent with the results of Chapter 4. It was determined that the mean height of the computational model was smaller than the Diffusion Model for chains with less than approximately 200 segments (see Figure 4.8). Dependence on the bending stiffness Similar simulations were performed for phantom and self-avoiding loops with constant length N = 90, but different values of the bending stiffness kc. A family of curves for the mean top height and gap distance are presented in Figure 5.7 for PPC100 loops. The SA100 data looked very similar, although the kc = 0 curve was not separated as strongly from the rest, likely due to the implicit bending stiffness resulting from sterics. 5. Single Loop under Parallel Confinement 88 Fig. 5.5: (a) Mean top height and (b) mean gap distance vs. compressive force for confined P P C 1 0 0 polymer loops of various lengths (N £ {15,30,45,60,75,90}). Force is scaled by kBT/b. The data is presented in a log-log plot. 5. Single Loop under Parallel Confinement 89 (a) /I I c 2.0 LU 1.5 n 1 r - - i 0.0 (b) £ S * 1 X • P P C 1 0 0 A S A 1 0 0 DM leading term ~ i • T" ~ i < n 20 40 60 80 100 0 20 40 60 80 100 Number of segments Number of segments Fig. 5.6: Power law prefactors and exponents of compressive force vs. number of segments for PPC100 and SA100 confined polymer loops. Various fits of F = a(h)b for the parameters a and b were performed for chains of different lengths in the strongly confined regime. kc — 0. Predictions from the leading term of the Diffusion Model (as the gap distance becomes small) are superimposed as a dashed line. 5. Single Loop under Parallel Confinement 90 1 0) 1 (b) TTT, | | I >—" H I I I 0.01 0.1 1 10 0.01 0.1 1 10 C o m p r e s s i v e f o r c e C o m p r e s s i v e f o r c e Fig. 5.7: (a) Mean top height and (b) mean gap distance vs. compressive force for confined PPC100 polymer loops of various bending stiffnesses (kc £ {0,1,3,7}). Force is scaled by k^T/b. The data is presented in a log-log plot. 5. Single Loop under Parallel Confinement 91 (a) • PPC100 * SA100 —1 ' 1 ' r 2 4 6 E f fec t i ve b e n d i n g s t i f fness coef f i c ien t (b) i • PPC100 * SA100 Ef fec t i ve b e n d i n g s t i f f ness coe f f i c ien t Fig. 5.8: Power law prefactors and exponents of compressive force vs. fceff for P P C 1 0 0 and SAfOO confined polymer loops. Various fits of F = a{h)b for the parameters a and b were performed for chains of different persis-tence lengths. N = 90. Power law fits were performed for force vs. gap distance in the strong confinement regime (F > 0.7). The results appear in Figure 5.8. Results are plotted against the effective bending stiffness keR, using the relationship determined in Chapter 3. Recall that within our system of units, fceff is simply equal to the persistence length for free chains in solution. The graph of the exponent illustrates the crossover from flexible chain to elastic rod-like behaviour. The exponent appeared to approach then stabilize to a value of 1.85(1) for fceff > 6 (this was the mean of all data points in this range). This was slightly larger than the value of 5/3 predicted for an untethered elastic rod. Consistent with theory, the prefactor also appeared to approach a stable value, loosing fceff dependence as the system switches from flexible to elastic rod-like.behaviour. This value was approximately equal to 55(3), as calculated by averaging the prefactors from the stiffest SA100 chain and PPC100 chain. It is also interesting to note that for fceff > 3, the phantom and self-avoiding data were roughly indistinguishable. This was consistent with previous observations that for a given contour length, as the persistence length increases, the relative size, and hence the importance, of the excluded 5. Single Loop under Parallel Confinement 92 (a) PPC100 SA100 — i 1 . 1 1 1 1 1 1 1 0 2 4 6 8 10 Effective bending stiffness coefficient (b) I PPC100 SA100 0 2 4 6 8 10 Effective bending stiffness coefficient Fig. 5.9: Power law prefactors and exponents of compressive force vs. fc0ff for P P C 1 0 0 and SA100 confined polymer loops. Various fits of F = a(h)b for the parameters a and b were performed for chains of different persis-tence lengths. N = 45. volume with respect to the effective segment length decreases. For comparison with the N = 90 data, the simulations were repeated for N = 45. The results of power law fits for F > 0.7 are presented in Figure 5.9. The graphs had the same characteristic shape as the N = 90 data. The prefactor and exponent approached the values of approximately 1.90(2) and 31(1), respectively (this was the mean of all data points with ^eff > 8). It was expected that the "elastic rod plateau" would be reached at approximately half the kes value required for the longer chain. While this appeared to be roughly true for the prefactor, it was less clear for the exponent. Further simulations on longer polymer chains could be performed to investigate this further. Finally, it is important to note that for the simulations with a non-zero explicit bending stiffness, the scaling exponents were different from the equivalent freely jointed chain. Thus, when studying compression un-der strong confinement, even for approximately 20 statistical segments, the freely jointed chain representation did not appear to be a good approxima-tion of the worm-like chain. 6. Concluding Discussion 93 6. C O N C L U D I N G DISCUSSION In this thesis, I have characterized how shape and elastic properties of a tethered polymer loop depend on the parameters of length, bending stiffness and excluded volume, plus parallel confinement. This was accomplished through the use of a Monte Carlo simulation. Properties in the large N limit were investigated and compared with analytical calculations. An explicit angular bending potential was included in the computational model, which allowed it to smoothly interpolate between the freely jointed chain (FJC) and worm-like chain (WLC) regimes. Rough bounds were placed on when a worm-like chain can be accurately represented by an equivalent freely jointed chain. In this concluding chapter, I will compare these findings with the results of a full network simulation performed by D. Boal [1]. Following this will be a discussion of the elastic thickness compressibility of a tethered polymer loop as found by simulations under parallel confinement. 6.1 Comparison with a Full Network Simulation The most notable simulations of tethered polymer networks were published as a computer model of the R B C spectrin network by David Boal [1]. In his simulation, each spectrin filament was represented by a chain of connected spheres, similar to the approach in this thesis. As here, steric interactions between the spheres were used to enforce self-avoidance in Boal's network. The filaments were connected at their ends in six-fold junctions to form a triangulated network. To represent the putative tethering with the lipid bilayer through band 3 proteins, Boal's network was "tacked down" to an impenetrable interface at the midpoint of each filament. Most important, Boal's network was placed in a periodic box with stress-free boundary con-ditions at the perimeter and top. Comparisons of my results for a single loop may be made with Boal's simulation by constructing an effective network representation. To do this, I used an affine approach, where the network was approximated by properties of the single loop in a triangulated array. For a network of equilateral trian-gles formed by connecting identical polymer loops end-to-end, the effective 6. Concluding Discussion 94 area per chain is simply Tge/2\/3, where vee is the distance between the end points of the polymer loop (also known as the end-to-end distance). Some notable differences between Boal's network simulation and the model used in the present study are the following: (i) In Boal's simulation, adjacent spheres were connected by string-like tethers (kt = 0 in the current nomen-clature), rather than by stiff spring-like tethers (kt = 100 in most of my simulations). However, results from Chapter 3 suggest that this does not have a large effect on shape properties, (ii) In my effective network model, the polymers are viewed as tethered at their end-points, rather than at their mid-points. While this may lead to different properties, it is not, in fact, clear from the literature which model best represents the architecture of the R B C . In this section I will use the effective network representation to derive three properties: the mass-average thickness, (z), the mean network surface density, (p), and the apparent thickness compressibility introduced by Boal for a transverse elastic modulus, Y [1]. A sketch of the R B C architecture and these geometric features are shown in Figure 1.1. The thickness z is a mass-weighted thickness of a chain configuration, whereas (z) represents the ensemble average of this quantity. To facilitate a direct comparison with the results in Boal's paper, this quantity is scaled by the contour length to form a dimensionless reciprocal thickness, L/(z). The network surface density is given by the number of chains per unit area. For the effective network model this is simply p = 2\/2>/r2e. The density is scaled by the density of the maximally stretched network, as defined by AC/(A), where (A) is the mean area per network junction, and Ac is the mean area/junction of the stretched out network. Y is an elastic constant that is based on the fluctuations of the mean thickness z. For the effective network model, it is defined as: Y = (p)kBT(z)/{(z2) - (z)2) (6.1) which becomes for a triangulated network, Y = 2V3kBT[(z)/((z2) - (z)2)]/(r2ee) (6.2) In Boal's simulation, this compressibility modulus was defined in a similar way, but instead scaled by the mean area of the periodic box rather than the mean area per chain. The corresponding dimensionless elastic constant was given as: (3YL3 (recall that (3 is simply l/kBT). Figure 6.1 shows the two measures of network geometry, L/(z) and AC/(A), as a function of the number of segments per chain. Within the range tested, the reciprocal thicknesses and surface densities of the PPC100 6. Concluding Discussion 95 (phantom) and SA100 (self-avoiding) single loop models predicted a thicker and more dense network than in Boal's simulation, fn other words, the single loop models were extended further from the interface and the chain ends were closer together on average. The reciprocal thickness of the SA100 simulation was roughly a factor of two smaller than in Boal's simulation. However, aside from this scaling factor, which appears in the log-log plot as a vertical offset, the behaviour of the SA100 and full network curves as a function of N was very similar. Both of these curves indicate a power law de-pendence on N above approximately 15 segments. The results of power law fits in this regime are presented in Table 6.1. The behaviour of the phantom PPC100 chain was very similar to the SA100 chain below 5 segments. For longer chains, the two curves rapidly diverged as the phantom chain obeyed a different power law for large N. Results of a power law fit for the PPC100 chain also appear in Table 6.1. The phantom chain did not grow away from the interface as quickly as the self-avoiding chains. The power law exponent for the reciprocal thickness was 0.417(8)1, which was slightly smaller than 1/2, the expected value for a long, phantom random flight chain. As inves-tigated in subsection 4.2.1 and Figure 4.8, approximately 200 segments were required to bring the PPC100 chain into the large N regime with respect to the measure of mean thickness. Model L/(z) Ac/(A) /Vi a b a b PPC100 2.15(6) 0.417(8) 1.4(1) 1.02(2) 15, 15 SA100 2.53(7) 0.289(8) 1.0(1) 0.83(3) 15, 15 Boal [1] 5.3 0.3 0.6 0.75 17, ? Tab. 6.1: Scaling laws for the reciprocal thickness L/(z) and the surface density (Ac/(A)) for tethered polymer networks. Da ta was fit to y = axb for N > A q . A plot of the surface densities for the three network models is presented in Figure 6.1 (b). Power law fits to each model are presented in Table 6.1. The results for surface density tell a complementary story to that of the thickness. As before, the SA100 effective network model and Boal's network model scaled very similarly with N, aside from a difference in the prefactor of roughly a factor of two. The density of the P P C 100 effective network grew at a faster rate than the self-avoiding network, with an exponent very close to the expected value of 1. For the self-avoiding networks models, it 1 The number is brackets indicates the uncertainty in the least significant digit 6. Concluding Discussion 96 1 10 100 1 10 100 Number of segments Number of segments Fig. 6.1: (a) Scaled reciprocal thickness (L/(z)) and (b) surface density (Ac/(A)) vs. number of segments (N) for tethered polymers in networks. The data is presented in a log-log plot. P P C 1 0 0 and SA100 refer to the effective net-work models that were based on single phantom and self-avoiding polymer loops. In these simulations there was no explicit bending stiffness (kc = 0) and stiff springs connected neighbouring spheres (kt = 100). The results of the Boal network simulation are overlaid [1]. A fit to the equation y = axb was performed on the P P C 1 0 0 and SA100 data for Ns > 15, and the results superimposed. The power law fit for the Boa l simulation as presented in reference [1] is superimposed as a dotted line. 6. Concluding Discussion 97 appeared that steric interactions within each chain caused it to swell. This increase in the end-to-end length likely caused a decrease in the mean height, commensurate with the outcome of thickness.2 From the Boal simulation, steric repulsions between neighbouring chains appeared to cause the network to swell even further in the lateral dimensions. However, the increased lateral swelling may have been caused primarily by the tethering of the filaments at their midpoints rather than their endpoints. Fluctuations in the mass-weighted thickness were used to compare com-pressibilities of these network models. First, a plot of the dimensionless material elasticity (3YL3 is given in Figure 6.2 along with results for the power law fits in Table 6.2. The error in this elastic constant was estimated at 10%. It was found that the effective network representation based on either PPC100 or SA100 chains was stiffer than Boal's network simulation. Within the range of N in Boal's simulations, the elastic constant for the net-work based on the SA100 chain was larger by a factor of two to three. This difference increased progressively with A^ due to the difference in power law exponents. Also, the network based on the PPC100 chain was even stiffer. Although the prefactor was approximately the same as for the SA100 model, the exponent for the PPC100 model was larger. The origin of these differ-ences appears for the most part to arise from the dependence on surface density in the three models. To see the effect of surface density, Figure 6.3 presents a comparison between elastic stiffnesses per chain for the single loop and Boal's network simulations. From this figure it is clear that the stiffness per chain was very similar for all three models. Lastly, it should be noted that a single elastic parameter like Y cannot be expected to capture the anharmonic response of soft tethered polymer loops. In order to show the true compliance of a tethered polymer loop, it was necessary to simulate the loop under parallel confinement over a large range in separation, as will be discussed in the next section. 6.2 Discussion on Elastic Thickness Compressibility The elastic thickness compressibility of a tethered polymer loop was exam-ined by simulations of a loop between parallel planes. A schematic of the confinement is given in Figure 1.1. The relationship between force, F, and gap distance, h, was studied, and the results compared to an analytical 2 In fact, the correlation between mean thickness and area is further explored in ref-erence [1], The volume per junction vertex was found to fluctuate very little for a fixed number of segments per chain. 6. Concluding Discussion 98 100-J 10-| , , — , . — . 1 10 100 Number of segments Fig. 6.2: Dimensionless elastic constant (/3YL3) vs. number of segments (N) for tethered polymer networks. Da ta is presented for the P P C f O O and SA100 effective network models, as well as for Boal network simulation [1]. In these simulations there was no explicit bending stiffness (kc = 0). A fit to the equation y = axb was performed on the P P C 1 0 0 and SA100 data for Ns > 15, and the results superimposed. The power law fit for the Boa l simulation as presented in reference [1] is superimposed as a dotted line. Model (3YL3 Ni a b PPC100 216(36) 1.58(5) 15 SA100 213(40) 1.23(6) 15 Boal [1] 146 1 ? Tab. 6.2: Scaling laws for the dimensionless elastic constant (3YL3 for tethered polymer networks. Da ta was fit to y = axb for N > N\. 6. Concluding Discussion 99 1000' co _J >-C O . CD JZ o 1 CD CL CO CO CD T3 CD CJ CO 100 H 10 • PPC100 * SA100 A Boal [1] - T 1—i—I—r— 10 —i—i—i—J— 100 Number of segments Fig. 6.3: Scaled stiffness per chain (3YLA/3 vs. number of segments (N) for teth-ered polymer networks. Data is presented for the PPCfOO and SA100 effective network models, as well as the Boal network simulation [1]. In these simulations there was no explicit bending stiffness (kc = 0). 6. Concluding Discussion 100 random flight model (the "Diffusion Model" (DM) described in Chapters 4 and 5). Figure 6.4 shows a comparison of the results for single self-avoiding (SA100) and phantom (PPC100) polymer loops, as well as the predictions of the Diffusion Model. Also superimposed are the ensemble averages of the mass-average height and the top height (i.e. the height of maximum excur-sion). The force response was highly non-linear, but exhibited the expected power law dependence under strong parallel confinement. Most difficult to examine, the steric repulsion at large gap distance was very weak, and was found to diminish precipitously above the mean top height of an unconfined loop. Under weak confinement, a harmonic spring was attached to the con-fining plate to act as a force transducer. Forces measured in this way showed good agreement with the random flight model when kc = 0. Here again, an equivalent network representation can be constructed from this single loop response to predict the dependence of normal stress on network thickness, which is shown in Figure 6.5. The stress is defined by the force divided by the area per chain, which is made dimensionless by a factor of kBT/b3. Under strong confinement, the force response was explored as a function of the number of segments, as well as bending stiffness, all of which was observed to obey power law scaling. It was found that the strongly confined behaviour of the PPC100 chain agreed with the predictions of the Diffusion Model, even for a chain with only 15 segments. When the bending stiffness was increased but the contour length held fixed, the polymer loop crossed over into the "stiff elastic rod" regime. This was marked by a change in the relationship between force and gap distance. Surprisingly, even when modelling a W L C with 20 statistical segments, the force response did not look like a freely jointed chain. This contrasts with the behaviour of most shape properties for the unconfined chain, where it was found that a W L C with 11 statistical segments could be well represented by an equivalent freely jointed chain. Thus, it may be important to include an explicit bending stiffness when studying the compressibility of a polymer network when each chain has less than 20 statistical segments. 6. Concluding Discussion 101 (a) P P C -m, 0.01 X 0.1 1 10 Compressive force (b) SA • Gap distance • Top height * Mean height DM gap distance nr, 10 Compressive force Fig. 6.4: Height measurements as a function of force for (a) PPCfOO and (b) SA100 confined polymer loops compared with the Diffusion Model (shown by dashed curves DM). Vertical error bars represent the standard deviation of the heights. Errors in the force are represented by horizontal error bars in the gap distance data. The loop was confined by two parallel planes with the untethered plane resting in a harmonic potential (HF simulations). 6. Concluding Discussion 102 (a) P P C 10 • Gap distance • Top height A Mean height TTTT, 0.01 rTTTl— 0.1 (b) SA Compressive stress 1E-3 0.01 0.1 Compressive stress Fig. 6.5: Height measurements as a function of stress for (a) PPC100 and (b) SA100 confined polymer loops. The stress was made dimensionless by a factor of ksT/b3. Vertical error bars represent the standard deviation of the heights. Errors in the stress are represented by horizontal error bars in the gap distance data. The loop was confined by two parallel planes with the untethered plane resting in a harmonic potential (HF simulations). Bibliography 103 Bibliography [1] D. H. Boal. Biophys. J., 67:521-529, 1994. [2] http://www.beyonddiscovery.org/beyond/beyonddiscovery.nsf/web/polymi [3] N. Mohandas and E. Evans. Annu. Rev. Biophys. Biomol. Struct., 23:787-818, 1994. [4] K. Svoboda, C. F. Schmidt, D. Branton, and S. M . Block. Biophys. J., 63:784-793, 1992. [5] B. T. Stokke, A. Mikkelsen, and A. Elgsaeter. Biophys. J., 49:319-327, 1986. As cited in [3]. [6] B. Alberts, D. Bray, J. Lewis, M. Raff, K. Roberts, and J. D. Watson. Molecular Biology of the Cell. Garland Publishing, Inc., New York, second edition, 1989. [7] V. Heinrich, K. Ritchie, N. Mohandas, and E. Evans. Biophys. J., 81:1452-1463, 2001. [8] E . Evans, K. Ritchie, and R. Merkel. Biophys. J., 68:2580-2587, 1995. [9] C. Picart and D. E. Discher. Actin protofilament orientation at the erythrocyte membrane. Biophys. J., 77:865-878, 1999. [10] P. G. de Gennes. Scaling Concepts in Polymer Physics. Cornell Uni-versity Press, Ithaca, New York, 1979. [11] P. Flory. Principles of Polymer Chemistry. Cornell University Press, Ithaca, New York, 1971. [12] J. D. Ferry. Viscoelastic Properties of Polymers. John Wiley and Sons, New York, second edition, 1970. [13] S. F. Edwards. Molecular Fluids, chapter The configurations and dy-namics of polymer chains. Gordon & Breach, London, 1976. [14] P. G. de Gennes. Phys. Lett, 38A:339, 1972. [15] J. des Cloiseaux. J. Phys., 36:281, 1975. [16] S. F. Edwards and K. F. Freed. J. Phys., A2:145-150, 1969. Bibliography 104 [17] A. K. Dolan and S. F. Edwards. Proc. R. Soc. Lond., A337:509-516, 1974. [18] M . Doi and S. F. Edwards. The Theory of Polymer Dynamics. Claren-don Press, Oxford, 1986. [19] F. C. MacKintosh, J. Kas, and P. A. Janmey. Phys. Rev. Lett, 75:4425-4428, 1995. [20] F. Gittes and F. C. MacKintosh. Phys. Rev., E58:1241, 1998. [21] E . Frey, K. Kroy, and J. Wilhelm. Advances in Structural Biology, vol-ume 5, chapter Viscoelasticity of Biopolymer Networks and Statistical Mechanics of Semiflexible Polymers, pages 135-168. JAI Press Inc., 1998. [22] E . Eisenriegler. Polymers Near Surfaces. World Scientific, Singapore, 1993. [23] D. H. Napper. Polymeric Stabilization of Colloidal Dispersions. Aca-demic Press, London, 1989. [24] D. H. Boal and U. Seifert. Phys. Rev. Lett, 69:3405-3408, 1992. [25] S. K. Boey, D. H. Boal, and D. E . Discher. Biophys. J., 75:1573-1583, 1998. [26] D. E . Discher, D. H. Boal, and S. K. Boey. Biophys. J. , 75:1584-1597, 1998. [27] D. Boal. Computational Physics. Simon Fraser University, Canada, 1997. [28] K. Binder and D. W. Heermann. Monte Carlo Simulation in Statistical Physics. Springer-Verlag, Berlin, 1992. [29] P. Stoltze. Introduction to simulation method, 1997. http://www.fysik.dtu.dk/ stoltze/text/simull/halle.html. [30] A. Lasota and M. C. Mackey. Chaos, fractals, and noise: Stochastic aspects of dynamics. Springer-Verlag, New York, second edition, 1994. [31] M . P. Allen and D. J. Tildesley. Computer Simulation of Liquids. Clarendon Press, Oxford, 1987. Bibliography 105 [32] G. R. Grimmett and D. R. Stirzaker. Probability and Random Processes. Clarendon Press, Oxford, 1992. [33] W. Rudin. Principles of Mathematical Analysis. McGraw-Hill, Inc., New York, 1976. [34] L. D. Landau and E. M . Lifshitz. Statistical Physics Part 1. Butterworth-Heinemann, Oxford, third edition, 1980. [35] E . Schrodinger. Statistical Thermodynamics. Dover Publications, Inc., New York, 1989. [36] P. J. Rossky, J. D. Doll, and H. L. Friedman. J. Chem. Phys., 69(10):4628-4633, 1978. [37] E. Evans and K. Ritchie. Biophys. J., 72(4):1541-1555, 1997. [38] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical Recipes in Fortran 90. Cambridge University Press, Cam-bridge, second edition, 1996. [39] P. J. Flory. Statistical Mechanics of Chain Molecules. John Wiley & Sons, Inc., New York, 1969. [40] O. Kratky and G. Porod. Rec. Trav. Chirn., 68:1106, 1949. As cited in [21]. [41] E . Frey, K. Kroy, and J. Wilhelm. The Wiley Polymer Networks Group Review Series, volume 2. 1998. arXiv:cond-mat/9808022. [42] W. Kuhn. Kolloid Z., 68(2), 1934. As cited in [18]. [43] P. J. Flory. J. Chem. Phys., 17:202, 1949. [44] P. E . Rouse. J. Chem. Phys, 21:1272, 1953. [45] W. Kuhn. Kolloid Z., 76:258, 1936. As cited in [39]. [46] W. Kuhn. Kolloid Z., 87:3, 1939. As cited in [39]. [47] P. Debye. J. Chem. Phys., 14:636, 1946. As cited in [39]. [48] H. Benoit and P. M . Doty. J. Chem. Phys., 87:958, 1953. As cited in [41]. [49] M . Fisher. J. Phys. Soc. Japan, 26 (Suppl.):44, 1969. As cited in [10]. Bi bliography 106 [50] K. Sole. J. Chem. Phys, 55:335, 1971. As cited in [53]. [51] J. Rudnick and G. Gaspari. J. Phys. A: Math. Gen, 19:L191, 1986. [52] J. A. Aronovitz and D. R. Nelson. J. Phys. France, 47:1445, 1986. [53] O. Jagodzinski, E . Eisenriegler, and K. Kremer. J. Phys. I France, 2:2243-2279, 1992. [54] L. D. Landau and E. M . Lifshitz. Mechanics. Butterworth-Heinemann, Oxford, third edition, 1976. [55] H. W. Diehl and E. Eisenriegler. J. Phys. A: Math. Gen., 22:L87, 1989. As cited in [53]. [56] G. Porod. Monatsh. Chem, 80:251, 1949. As cited in [39]. [57] R. Dennemeyer. Introduction to partial differential equations and boundary value problems. McGraw-Hill Book Company, New York, 1968. [58] H. S. Carslaw and J. C. Jaeger. Conduction of heat in solids. Clarendon Press, Oxford, 1959. [59] Q. McNemar. Psychological statistics. John Wiley and Sons, Inc., New York, fourth edition, 1969. [60] E . W. Weisstein. The CRC Concise Encyclopedia of Mathematics. CRC Press, USA, 1999. [61] P. M. Morse and H. Feshbach. Methods of Theoretical Physics. McGraw-Hill, New York, 1953. [62] E. Evans, 2000. personal communication. Appendices 107 APPENDICES A. Glossary of Symbols 108 A. GLOSSARY OF SYMBOLS Number of vertices N Number of segments (Ny — 1) b Bond (or segment) length1 L Contour length (N x b) Radius of the spheres h Intervertex spring coefficient kc Bending stiffness coefficient fceff Effective bending stiffness F Compressive force DT Step parameter for vertices DTB Step parameter for the upper plane fcB Boltzmann constant T Temperature l/kBT (A) Ensemble average of the quantity A ree End-to-end distance r9 Radius of gyration M Mean <7 Standard deviation Z Mass-average thickness t Top height ML Top height elastic modulus Y Pseudo transverse Young's modulus h Gap distance 1 Occasional ly , b signifies the exponent of a power law fit. It wil l be clear f rom context when this is the case. A. Glossary of Symbols 109 r.v. Random variable R B C Red blood cell SMC Smart Monte Carlo M C Standard Monte Carlo F J C Freely jointed chain W L C Worm-like chain SA[/V] Self-avoiding Chain computational model [with kt = N] P P C [TV] Partial Phantom Chain computational model [with kt = N] D M Diffusion Model (random flight model) SA SS Self-avoiding "smart scaling" HF Harmonic compressive force C F Constant compressive force B. Basic Concepts of the Monte Carlo Method 110 B. BASIC C O N C E P T S O F T H E M O N T E C A R L O M E T H O D This appendix discusses the basic ideas and theory behind the Monte Carlo method. It is essentially an expansion of the material in Section 2.1. The most significant addition is a section which summarizes the relevant formal-ism of the mathematical theory of probability. However, specific implemen-tations of the Monte Carlo method, and its application to physics are not addressed here. For a more practical mini-tutorial of this computational technique, the appendix should be read together with Section 2.2. The Monte Carlo Method refers to a collection of stochastic techniques that were developed by Metropolis, Ulam, and von Neumann in the 1940s at Los Alamos. For this reason it is sometimes called the Metropolis method. One of the central problems that this method addressed was the numerical calculation of integrals through the use of statistical sampling. This is the technique that will be presented here. The discussion is based primarily on references [27], [28] and [29]. Consider the problem of approximating the 1-dimensional integral: If this integral is well-defined, it can be approximated in the standard way by selecting N equally spaced points in the interval [0,d], then evaluating the function / at these locations. This yields the expression, B.l Simple sampling (B.l) L E , " i / ( * i ) A (B.2) app — N where S = {x^} is the set of uniform x coordinates, and A is the normaliza-tion factor d. Clearly, as N approaches oo, Iapp approaches / . B. Basic Concepts of the Monte Carlo Method 111 The idea behind the Metropolis approach is this: the x^s do not have to be chosen in a grid in order to yield a meaningful approximation; rather, they can be chosen randomly between 0 and d. Then if the new set [xi] is selected with uniform probability within this interval1, equation (B.2) still represents a valid approximation of the integral / as N becomes large. Now clearly, this technique is not limited to integration over one dimen-sional intervals, but can be easily extended to more general domains. For if we call the domain of integration fi, and choose the x^s uniformly in fi, then equation (B.2) remains our approximation. However, the normalization factor is no longer simply A = d, but is represented more generally by the integral over all of fi (in simple cases, the "volume" of fi). The domain of integration can be quite general, which makes the Monte Carlo integration method both versatile and powerful. In many cases, this approach more appropriate than discretizing fi into a grid. In order to state these ideas precisely, I will present a little bit of formalism in the next section. B.2 Formalism: Probability spaces, random variables, and distribution functions The discussion may be clarified by introducing some of the mathematics behind the Monte Carlo method. The definitions are based primarily on reference [32], G. R. Grimmett and D. R. Stirzaker. However, they will occasionally be more intuitive than mathematically complete, and the in-terested reader may wish to peruse the references at chapter's end. On the other hand, readers who are already familiar with basic probability theory may wish to skip this section. B.2.1 Probability spaces Consider a random experiment whose outcome is unpredictable. We call the set of all possible outcomes the sample space, denoted by the symbol fi. Subsets of the sample space which form a mathematical structure called a sigma field are called events.2 Events may be thought of as "things that can occur when the experiment is performed". There is a probability asso-ciated with the occurrence of each event, denoted by the symbol P(event). The function which tells us the probabilities of different events is called the 1 This will be clarified in the next section. 2 Specifically, a sigma field T (also known as a sigma algebra) is a collection of subsets of Q such that: i) Q G T, ii) / € T implies that fc G where fc is the complement of / , and iii) (JXi /» ^ 3~< f ° r a " /» ^ -^i *"n a^ ' s ' s dosed under countable unions. B. Basic Concepts of the Monte Carlo Method 112 probability measure. Let A be an event. Then according to the standard interpretation of probability, P(A) simply represents the fraction of times that A will occur if the experiment is repeated many times. Taken together, these three fundamental elements—the sample space (fi), the set of events (the sigma field, J7), and the probability measure (P)—form a structure that is called a Probability Space, W = (fi,.?7, P). Let A and B be events in the probability space W (i.e. they are subsets of fi which are elements of IF). Then the conditional probability that A occurs given that B occurs is given by Now consider a function that maps from the sample space to the real num-bers, e.g. X : fi —> K. If for all x, the set {LO € fi : X(to) < x} is an element of the sigma field J7, then X is called a random variable on the probability space W. X takes on different values with a probability determined by the measure P. It can be thought of as a measurement that is made on the outcome of the experiment. Specific realizations of X, or "individual mea-surements", are called instances of X. Making such measurements to create a set of instances is called sampling. Now in a more general context, X need not map to the real line, but may map to any set with a sigma algebra. This set is called the state space of X. Thus, measurements can be more than simply numbers. However, to simplify the current discussion I will focus on random variables (r.v.'s) which map to the reals. Rather than always dealing directly with the probability space, a ran-dom variable can be characterized by its distribution function, or cumulative distribution. The cumulative distribution, Dx : K —> [0,1], gives the prob-ability that the outcome of the experiment will be mapped to a value less than or equal to x by the function X: B.2.2 Conditional probability P{A\B) = P(AC)B) P{B) (B.3) where A n B is the intersection of sets A and B. B.2.3 Random variables and distribution functions (B.4) Here, the event {x < X} is understood to mean the set of events {LO : X(u>) < x}. Although it is understood that X simply maps random events to the real B. Basic Concepts of the Monte Carlo Method 113 numbers, I will occasionally adopt the common attitude of "downplaying the probability space", by speaking as if X itself takes on random values. Now if Dx{x) is differentiable, its derivative is called the probability density of X, and denote it by p(x). It is clear how it gets its name, for p(x)dx represents the probability that X maps to values in the interval \x,x + dx]. If p(x) is constant over some range of R, then X is said to be uniform or to have a uniform distribution in this range. Now clearly, this discussion of p(x) does not apply to discrete random variables, which do not map from a continuous state space. Such variables are generally characterized by simple jump discontinuities in Dx(x). In these cases, the role of p(x) is replaced by the probability mass function, p(x), which gives the probability that X will take the value x (i.e. p(x) — P({u) £ fi : X(oj) = x})). However, this can be accommodated this within the probability density formalism by representing the derivative of a jump discontinuities of height by uk x 5(x — xn)" • Here, 5(x) represents the standard delta function, and XQ is the position of the step. In any case, it is clear that P(x e (a, b}) = (B.5) These definitions extend quite naturally to multiple random variables. Consider the random vector, X : fi —> R N , consisting of the N random variables, Xi : fi —> IR. Then the distribution function becomes Tx • ^ N —> [0,1], with Tx(x) = P ( X 1 < x 1 , . . . , X N < x N ) (B.6) where X j is the ith component of x. Again, if Tx is differentiable, the proba-bility density is given by the Jacobian of Tx (determinant of the derivative), and is denoted p{x\,..., x^) as before. Analogously, p(x\,... ,x^)dx\... dxjy represents the probability of X mapping to the cube [xi,xi + dx\] x . . . x [xN,xN +dxN}. When the sample space is fully characterized by N random variables, the old probability space can be replaced with a new one where fi = RN, T is the standard Borel sigma field, and the measure is induced by the joint probability of the random variables, Tx- In this case, we can apply the discussion above, and talk about the probability density, and distribution function of the sample space itself. B. Basic Concepts of the Monte Carlo Method 114 B.2.4 Expectation values and functional composition The expectation value of a random variable is also known as the mean, or average. In the standard way, it is given by: wen This can be written in a way which also looks good for continuous variables, by using measure theory notation (see reference [33] for more details): (X) = / X(w)dP(w) (B.8) Jn However, in either notation, the meaning of the average is clear. Now if the distribution function of a random variable is known, Tx, this can be used instead to define the mean: /+ O0 xdTx(x) (B.9) -oo where Tx is said to induce a measure on R (or to act as a Stieltjes function). If the density function exists, we can write: /+oo xp{x)dx (B.IO) -oo where the integration is done with the standard Lebesgue measure (normal integration). Given a function / : R —> R , we can compose it with X to construct a new random variable: f(X) : Q, —> R . Then, the average value of f(X) can be calculated as above, e.g. /+ O0 f(x)p(x)dx (B . l l ) -oo B.2.5 The law of large numbers In a practical sense, an experiment can only be performed a finite number of times, and so the true distribution function of a measurement (or random variable) cannot be determined. Fortunately, however, there are theorems which state that as the sample size increases, the difference between the measured distribution and the "true" distribution decrease. I will not go B. Basic Concepts of the Monte Carlo Method 115 into the details of such theorems here, but state the practical consequence that will be used. Consider [X] = (B.12) where { X J } consists of N instances of X. Then according to the law of large numbers, as TV approaches infinity, the error \[X] — {X)\ approaches zero. Though there are cases where this will not work, the truth of this theorem is broad enough for the current study. See reference [32] for further information. B.3 Importance Sampling So the essence of the M C integration method is to interpret an integral as the average value of some random variable, and then to approximate this average by sampling the random variable, as in equation (B.12). As an example, consider again the one-dimensional integral given in equation (B.l) . Let X be a random variable with probability density p(x). Furthermore, let us require that p(x) is non-zero only on the domain of integration, [0, d\. Next, compose X with the function: n(x) = 4^ r, for x s.t. p(x) ^ 0 = 0 otherwise (B.13) The average value of the r.v. g{X) yields the desired integral: /+oo g(x)p(x)dx -co Jo P{x) d f(x)dx (B.14) Then we can approximate this integral by averaging specific instances of g{X), with the probability distribution of X given by p(x). Notice that dividing by p takes care of the normalization. When the distribution of X is not uniform over the interval [0,d], i.e. p(x) is not a constant over this domain, the technique is called importance sampling. This method can be more efficient than "simple sampling" as p B. Basic Concepts of the Monte Carlo Method 116 can be chosen to preferentially sample regions of the domain of integration where f(x) is the greatest. So to summarize the basic strategy: to integrate a function / over some domain, we pick a random variable that is non-zero only on this domain. Then the integral is the average value of (f /p)(X), where p is the probability density of X. This average can be approximated by sampling the random variable. B.4 The Markov Chain In order to carry out the Monte Carlo method, it is necessary to construct a random variable with a specified probability density. A powerful technique for accomplishing this was developed by Metropolis et al. The basic idea is to form a sequence of random variables, specifically a Markov Chain, which converges to the desired r.v. In this section, I will present a brief introduction to Markov Chains, followed by a description of the prescription. Consider a sequence of random variables indexed by a "time variable": {Xi : i £ Z>o}. To simplify the present discussion, it will be assumed that the index is discrete. Furthermore, let us also assume that the random variables are discrete, and map to the (ordered) state space S. Let the probability mass function of each variable Xi be denoted by pi. Now, if each r.v. depends on the previous variables in a well-defined manner this sequence is called a random process. If in addition, the (conditional) distribution of Xi depends only on Xi-\, the random process is called a Markov Chain. Formally, this requirement is P{Xn = S\XQ = xo,..., Xn_\ — XN-\) = P(Xn — s\Xn-\ = xn-\) Vn > 1 and all s, XQ, . . . , x n - \ e S (B.15) Now suppose we want to construct a random variable X with probability mass function p. The prescription given by Metropolis et al. yields a Markov Chain such that for sufficiently large A'', X^ has a distribution which is sufficiently close to p. The key to this construction is the condition of detailed balance: p(i)Wi^f=p(f)Wf_i, V t , / G S (B.16) where Wa^b represents the transition probability from state a to state b, that is Wa^b = P(Xk = b\Xk^ =o)Vfc (B.17) When this condition is satisfied, we would like to have the following conver-gence: for all e > 0, there exists an N > 0 such that | | -P(X n = x) — P(X = B. Basic Concepts of the Monte Carlo Method 117 x)\\ < e, whenever n > N. Although there are some subtle issues which will not be addressed here, for practical purposes this convergence can usually be taken to hold true. A simple argument will now be presented to show that this is reasonable. B.4.1 Heuristic argument Here is a fairly well-known heuristic argument that a Markov Chain which satisfies condition (B.16) has the desired properties. It is based on a discus-sion in reference [28]. Consider a very large collection of M Markov Chains in parallel. We may think of this as a single Markov Chain, where each random variable in the sequence is now a random vector of M randon variables. Now suppose that the state space, S, of each random variable is simply {i,f}. Let the number of occurrences of the states i and / in an instance of the kth random vector, be given by ik and fk, respectively. Since M is large, the change in %k from the fcth instance to the (fc + l)th instance is given by: 5i = i k + i - ik « Wf->ifk — Wj_>/ife (law of large numbers) = ( B 1 8 ) Furthermore, by the condition of detailed balance, we have: *-*xW'A£-m) (ai9) Now ik x Wf-+i is a positive number, so it is clear that 5i is positive if the ratio fk/ik is greater than the desired ratio p(f)/p{i), and negative if fk/h < p{f)/p(i)- Thus, if the state i is under-represented in the instance of the kth random vector, the number of occurrences of i increases on average in the next instance of the Markov Chain. On the other hand, if i is over-represented, the number of i occurrences will decrease on average. The same holds true for /fc. So it appears that the Markov Chain approaches the desired distribution. Now this argument is not very satisfying for a number of reasons. Pri-marily, it only considers state spaces of two states, which are rarely of in-terest. Secondly, it does not consider the possibility of cycles, or other such phenomena. However, this argument provides a way to think about the problem, and can help us to develop some intuition about the Markov Chain. B. Basic Concepts of the Monte Carlo Method B.4.2 Algorithm In practice, one can generate an instance of the desired Markov Chain in the following way: 1. Initialization. Choose a starting state XQ = a, where Xi refers to the instance of the r.v. AV Also, initialize the index: i = 0. 2. Make a trial move. Consider the trial move as an instance of the ran-dom variable X t r i a i - The conditional probability of choosing Xtrial = / given Xi = xt is denoted by P ( X t r i a l = / | X c u r r e n t = x{). We write ^current instead of Xi to emphasize the fact that this probability is usually chosen so that it does not depend explicitly on i 3. Choose or reject this move according to some prescribed probability, P(Accept move? | X c u r r e n t = Xi & Xtrial = / ) • 3 As before, we are assuming that this probability does not depend explicitly on the index i. If the move is accepted, let X J + I = / . Otherwise, let Xi+\ = Xi. 4. Increase i by one and go back to step 2. Now this procedure gives us a lot of freedom both in choosing a trial move, and in setting the probability of acceptance. The primary condition that must be satisfied is detailed balance (equation B.16). In this context, Wa^b = P ( X t r i a i = 6 | X c u r r e n t = a) x (B.20) P (Accept move? | X c u r r e n t = a & A" t r i a i = b) Now the usual practice is to choose the trial moves symmetrically, so that P(^tr ia l = / ^ c u r r e n t = i) = P(XtTia\ = i |A" c u r r e n t = / ) • Then detailed balance can be satisfied by choosing the acceptance probabilities appropri-ately. One standard choice is accept a move from state i to state / with probability: minimum(l,p(/)/p(i)). Here, p(a) refers to the probability of finding the system in state a. An alternative choice for these probabilities is presented in section 2.2.4. Another important requirement, is that the probabilities be chosen so that the Markov Chain can access all states in S from any starting state. If these two conditions hold, then for sufficiently large N, XN should be essentially independent of the starting configuration x$. Thus, recording the 3 To be explicit, P(Accept move?) is shorthand for P(Accept move? = yes) where "Accept move?" is a random variable with state state space {yes, no}. Similarly, P(Reject move?) is equivalent to P(Accept move? = no) B. Basic Concepts of the Monte Carlo Method 119 state of the Markov Chain after every TV steps is equivalent to sampling the random variable XN repeatedly 4 . This is what is usually done in practice— sample a single Markov Chain every N steps. B.4.3 Continuous random variables Up to this point, I have discussed Markov Chains in the context of discrete random variables. This has served to make the discussion and theory rela-tively simple. However, these ideas can be extended to continuous random variables in a natural way. This topic will be discussed briefly in the current section, but continuous time/indexing variables and pathological functions will not be addressed. Suppose we wish to construct a random variable X : fi —> R, with probability density p(x). As before, the critical condition for generating a Markov Chain which converges to X is detailed balance. However in this context, the relevant states (subsets of the state space) are intervals (1-d volumes). That is, the condition of detailed balance can be written as: W i = [ a A _ f = { c 4 ] - (B.21) Jc p{x)dx Now in practice, one can usually carry out the procedure described above for discrete random variables by simply replacing the probability mass func-tions p with the probability density functions p. This corresponds to treating the state X = a 6 R, as the infinitesimal interval [a,a + dx]. Care must be taken, however, when using trial moves which map between unequal vol-umes of the state space. One common example of this is a scaling move, in which the current state is dilated or contracted. For the sake of concrete-ness, suppose that the initial state is i = x, and the trial move / is scaling by g. Then for detailed balance to be satisfied, we must have: ra+e rga+ge / p{x)dxW^f= / p(x)dxWf^i (B.22) J a J ga which implies that p(a)W^f = gp{ga)Wf^i (B.23) Notice that one' must include a factor of g when calculating the probability of accepting or rejecting this move. An application of this idea is discussed on page 124 of reference [31]. 4 In fact, the desired condition is that the process is ergodic. For more on ergodicity see reference [30]
Thesis/Dissertation
eng
Physics
Vancouver : University of British Columbia Library
For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use.
Graduate
Modelling a tethered polymer loop : shape and thickness compressibility
Text