Conformation-dependent sequence design of polymer chains in melts
Elena N. Govorun, Ruslan M. Shupanov, Sophia A. Pavlenko, Alexei R. Khokhlov
11 Conformation-dependent sequence design of polymer chains in melts
Elena N Govorun,* Ruslan M Shupanov, Sophia A Pavlenko, Alexei R Khokhlov
Faculty of Physics, M. V. Lomonosov Moscow State University, Leninskie gory, Moscow 119991 Russia E-mail: [email protected]
Abstract
Conformation-dependent design of polymer sequences can be considered as a tool to control macromolecular self-assembly. We consider the monomer unit sequences created via the modification of polymers in a homogeneous melt in accordance with the spatial positions of the monomer units. The geometrical patterns of lamellae, hexagonally packed cylinders, and balls arranged in a body-centered cubic lattice are considered as typical microphase-separated morphologies of block copolymers. Random trajectories of polymer chains are described by the diffusion-type equations and, in parallel, simulated in the computer modeling. The probability distributions of block length k , which are analogous to the first-passage probabilities, are calculated analytically and determined from the computer simulations. In any domain, the probability distribution can be described by the asymptote ~ k at moderate values of k if the spatial size of the block is less than the smallest characteristic size of the domain. For large blocks, the exponential asymptote exp( const dka ) is valid, d as being the asymptotic domain length ( a is the monomer unit size). The number average block lengths and their dispersities change linearly with the block length for lamellae, cylinders, and balls, when the domain is characterized by a single characteristic size. If the domain is described by more than one size, the number average block length can grow nonlinearly with the domain sizes and the length d as can depend on all of them. Keywords: polymer sequence design, first-passage probability, random multiblock copolymers, computer simulations, diffusion equations
1. Introduction
Polymer sequences can control the conformation of single macromolecules, as well as the spatial structure of polymer solutions and melts. This ability is obviously revealed for biomacromolecules with the unique native structure of proteins and the specific conformations of DNA and RNA [1-4]. The attempts to design artificial biomimetic polymers have led, in particular, to the idea of the conformation-dependent design of protein-like copolymers with hydrophobic core and polar shell. In a homopolymer globule consisting of hydrophobic monomer units, the monomer units in the surface layer are transformed into polar ones. The macromolecules thus obtained possess the ability to form soluble globules, which are stable against aggregation with each other. First protein-like copolymers were designed and described in the computer experiment using the Monte Carlo method [5]. Then, the statistical properties of sequences and the globule stability were investigated in the theory and computer simulations [6-8], and the higher stability of polymer globules formed after the conformation-dependent chemical modification was demonstrated experimentally [9, 10]. Besides, the evolutionary conformation-dependent process of mutations (monomer unit exchange) was studied in the theory and computer simulations, the final macromolecular conformations were hydrophobic-polar core-shell structures and core-tail (tadpole) conformations [11]. In addition to the sequences of core-shell globules, the conformation-dependent sequence design was theoretically described for the macromolecules adsorbed at a penetrable liquid-liquid interface [12]. For homopolymers adsorbed at a solid plane surface, the conformation-dependent sequence modification was carried out in the computer simulations, and the adsorption of macromolecules with the different types of monomer unit distributions was studied both in silico and experimentally [13, 14]. Aside from single macromolecules, all chains of a melted polymer can be modified. It is microphase separation that plays an important role in the material properties of block copolymers. The simplest lamellar structure was used as a pattern for the conformation-dependent sequence design in the homopolymer melt, and the statistical characteristics of sequences were studied in the theory and computer simulations. Using the dissipative particle dynamics, the lamellar-type structure was reproduced in the melt of the designed copolymers after switching on the repulsive interactions between the monomer units of different types [15]. The macromolecules obtained via the conformation-dependent modification can be considered as random multiblock copolymers, where a block is a sequence of monomer units of the same type bounded by the monomer units of another type. The monomer unit sequences can be characterized by the average composition and the block length distributions. At present, a great progress in the synthetic methods for multiblock copolymers with controlled length and degree of blockiness attracts attention to these materials [16-20]. In a large homogeneous globule and in a melt, the polymer chain conformations are random, and the dependence of the monomer coordinates on its number along the chain can be reduced to the trajectory of a random walker. The statistical weight of chain trajectories can be described by the Green function obeying the diffusion-type equation in terms of the Gaussian chain model [21]. For hydrophobic blocks in the globule core of a protein-like copolymer, the block termination at the core surface is determined by the absorbing boundary conditions for the diffusion equation [6]. For copolymer blocks corresponding to the lamellar pattern in melts, the block length distribution is also determined by the absorbing boundary conditions. In the probability theory, such analysis relates to the classical first-passage problem [22, 23]. Generalization of the problem to the cases of the moving walls or the random walks of special types continues to attract the interest [24-29]. In the present paper, the block length distribution and the number average block length correspond to the first-passage probability and the first-passage time in the 3D domains of different shapes. For the polymer blocks in a globule core or in a lamella [6, 15], the asymptote for the distribution of block length, k , is ~ k at moderate values of k , and such a power-law distribution can be described in terms of the Levy-flight statistics [30]. At very large values of k , the probability distribution of block lengths decreases exponentially. For protein-like copolymers, the wide block length distributions lead to the existence of long-range correlations in the sequences [6]. In the present paper, we put together the descriptions of sequence statistics for the alternating lamellae and periodically arranged cylinders and balls. The most detailed analysis is performed for the first considered pattern of hexagonally packed cylinders. In the section Model and methods , the system, the patterns, and the model approximations are described, the diffusion equation with the boundary conditions and the asymptotes are written out, and the dissipative particle dynamics (DPD) model is presented. In the section
Results and Discussion , the exact solutions of the diffusion equation in the considered domains are written out, the smoothing procedure for the computer simulation data is explained, the characteristics of the block length distributions are calculated, and their peculiarities are analyzed. In the Conclusions section, the results of the work are summarized.
2. Model and methods
Initially, we consider a homogeneous melt of homopolymer molecules each consisting of N monomer units. The polymer chains are assumed to be ideal and their conformations are described as random walks with the mean-square distance a between two neighbor monomer units in a chain. In that case, the distribution of distances between any two monomer units of the same polymer chain is described by the Gaussian statistics [21]. Then, we introduce an imaginary two-component bulk pattern that resembles one of the typical morphologies of microphase separation in block copolymer melts. Those morphologies are alternating lamellae, hexagonally arranged cylinders, and balls forming a cubic lattice, in particular, a body-centered cubic lattice (figure 1). Type A is assigned to the monomer units in the yellow (light) lamellae, cylinders and balls, whereas type B is assigned to all other monomer units in the blue (dark) domains. The sequence of morphologies in figure 1 follows a decrease in the volume fraction of yellow domains and, correspondingly, decreasing in the mean fraction of monomer A units in a polymer chain. When moving along a polymer chain, the type of monomer units changes at the interfaces between blue and yellow domains. Thus designed chains can be described as random multiblock copolymers, the chains of which consist of blocks, i.e., sequences of chemically identical monomer units. Figure 1. Scheme of modification of an initially homopolymer melt into a patterned melt of multiblock copolymer. Monomer units in the yellow (light) lamellae, in the cylinders and spheres are of type A, other monomer units are of type B. The cross-sections of lamellae (a), hexagonally packed cylinders (b), and a body-centered cubic lattice of balls (c) are shown. For the last pattern (c), the projections of the inner balls of elementary cells onto the section plane are shown by hatching. Monomer unit distribution along a polymer chain can be completely described by the distributions of A and B blocks over their length together with the average chain composition. In the work, those statistical characteristics of monomer unit sequences are determined via two independent methods. The first one is based on the calculation of the Green function of a polymer chain in the considered domains via solving the diffusion-type equations. The second one uses computer simulations of a homogeneous polymer melt with the geometrical patterns (figure 1), where the number of blocks of each length is directly counted.
The use of the diffusion-type equations for the description of polymer conformational statistics is based on the analogy between a random walk trajectory and a dependence of the position of a monomer unit on its curvilinear coordinate along an ideal chain. The statistical weight of chain conformations with fixed ends is described by the Green function of a polymer chain [21]. Polymer block trajectories begin at interfaces between A (yellow) and B (blue) domains and finish somewhere at those interfaces. Further, a domain (or monomer unit) type is denoted by ; = A,B. The Green function G ( r , k | r ) of the sequences containing k monomer units of type , beginning at the point r and ending at the point r in the domain obeys the equation )|,(6)|,( rrrr kGakGk , V , rr , (1) where is the Laplacian with respect to r . The absorbing boundary condition describes the sequence termination if its trajectory hits the surface S of domain: S kG r rr . (2) The initial condition is )()|0,( rrrr G (3) A part of the polymer chain in the domain with both ends at the interfaces becomes an block, and then the initial point r should be taken near the interface. Note that macromolecules are assumed to be long enough, so that the probability of chain termination inside a domain can be neglected. This condition is hold for Na >> d , = A,B , where d is the characteristic size of domain. For sequences of length k with the characteristic spatial size ak being much smaller than a characteristic domain size d , the Green function can be approximated by the solution of Eq. (1) in a semi-space x ),|,()|,()|,()|,( zkzGykyGxkxGkG rr (4) where |, ykyG )(2 3exp2 3 yykaka is a univariate Gaussian distribution and )|,( xkxG )|,()|,( xkxGxkxG . In this approximation, the Green function (4) does not depend on the segment type . The characteristic domain size d = L for lamellae, d A = R for cylinders, and d A = R s and for balls (figure1). We are aimed to determine the probability distributions P A ( k ) and P B ( k ) of A and B blocks, respectively, over their length k . The product P ( k ) k should give the fraction of blocks with the length in the interval [ k , k + k ]. It can be related to the integral V kGd |, rrr , which is equal to the number of chain trajectories consisting of k steps, beginning at the point r and staying inside the domain, in proportion to the total number of possible trajectories. Then, the probability distribution P ( k ) is equal to V kGdkkP |,)( rrr (5) or, with the use of Eq. (1), it can be calculated as an integral probability flux through the surface: S GdsakP n , (6) where the unit vector n is perpendicular to the surface. The absorbing boundary condition leads to G
0 at k . Then, from the formula (5) and the initial condition (3), the block size distribution P ( k ) should obey the normalization condition: kPdk A (7) For the solution (4) in a semi-space, the distribution of blocks over their length can be found with the use of Eq. (6): kaxkxkax xkxGakP x (8) This dependence can be interpreted also as a survival function describing the time of first passage of the plane x = 0 by the random walk [22]. A sequence of blocks in a polymer chain can be considered as a unilateral walk with the step size distribution described by the alternating functions P A ( k ) and P B ( k ). For the value x being of the order of a and for k >> 1, the dependence of P ( k ) on k (Eq. (8)) can be approximated as P ( k ) ~ 1/ k . 1 « k « ( d / a ) (9) A random walk with such a power-low distribution of step size (Eq. (9)) is often referred to as a Levy flight [30]. For larger blocks of the spatial size that is comparable with the characteristic size of the domain d , this size and the shape of the interface become substantial. For very long blocks with trajectories avoiding the boundaries, both the Green function and the block length distribution decrease exponentially with k . Having in mind that the distribution P ( k ) should be continuously extended over k > k * = ( d / a ) , the asymptotic form exp~)( dkadakP , k » ( d / a ) (10) can be written out, where the coefficient depends on the geometrical parameters of the region, the value of being of the order of unity. It is worth noting that the asymptotic form (10) describes the distribution P ( k ) for the case of a single characteristic domain size d . If the shape of the domain is complex and it is characterized by several spatial sizes, the behavior of the function P ( k ) at large k could depend on several parameters. The power-low dependence (9) at small k is valid until the block spatial size ak becomes comparable with the smallest characteristic size of the domain. Useful characteristics of a block length distribution are the number average length n k , mass average length w k , and dispersity Ð : kkPdkk n , kPkdkkk nw , Ð nw kk (11) The number average length is analogous to the first-passage time of a random walker. For the asymptotic forms (9) and (10), the dependences (11) of the average lengths on the characteristic size d of the domain are .~,~ adkadk wn (12) Then, the dispersity of the block length distributions Ð ad ~ (13) grows linearly with the domain size. The exact solutions of Eq. (1) are known for a ball, a cylinder, and a lamella [30], whereas such solutions are unknown for an infinite domain with the excepted domains in the shape of cylinders or spheres (as the blue domains in figure 1 (b) and (c)). In diblock copolymer melts, the sizes of the domains of different types are usually comparable. In this case, the probability of trajectories that never hit absorbing boundaries is negligible. In the theoretical analysis, we approximate these infinite domains by the simplest finite ones. For B blocks beginning near a certain cylinder (the pattern (b) in figure 1), we approximate the domain by a cylindrical layer surrounding this cylinder. For B blocks with the beginning near a certain ball (the pattern (c)), the approximation of a spherical layer is used. For some of the domain shapes in figure 1, the block length distributions for the pattern-modified polymers were calculated previously. Blocks in a ball can be described by the distribution function derived for hydrophobic blocks in a spherical core of a large globule formed by the protein-like copolymer [6]. Besides, the distribution of the polar block lengths in a protein-like copolymer was calculated [6], where an additional reflecting boundary condition at the outer surface was used for the spherical layer around the core, not as in the present consideration for B blocks. The block length distribution for a lamellar structure was written out and modeled in the computer simulations in the work [15]. In computer simulations, we generate a homogeneous polymer melt, assign types A and B to the monomer units according to the geometrical patterns in figure 1, and find the block length distributions. Polymer chains are modeled as identical phantom beads (monomer units) connected by “springs”. Initially, polymer chains are randomly distributed in the simulation box with periodic boundary conditions and the conformations of polymer chains are set as random walks with a constant spatial step. Then, the system of beads is equilibrated using the method of dissipative particle dynamics [32] until the conformational characteristics of polymer chains stop changing. The coordinates of beads are described by the radius-vectors { r i }. The force f i acting on the i -th bead in the system is a sum of the bond stretching forces b ij f acting from the neighbor beads in the chain, the conservative forces c ij f , dissipative forces d ij f , and random forces r ij f : ij ijijijiji rdcb fffff (14) In a Gaussian chain, the bond stretching force depends linearly on the distance between beads i and j : ijij K rf b , jiij rrr (15) where the constant K characterizes the bond stiffness. The forces c ij f , d ij f , and r ij f are acting from all j -th beads in the sphere of the cut-off radius r c with the i -bead in the center. The conservative forces c ij f are taken in the form c ijijijijijij rrrr rf (16) where is the interaction parameter. The bead mass m , the cut-off radius r c , and the thermal energy k B T are taken as reference units. Expressed in those units, the stiffness coefficient K = 4, the initial bond length is equal to 1, the interaction coefficient = 25, and the mean number density of monomer units is equal to 3. The dissipative forces d ij f and random forces r ij f are set as in the work [33] to provide the value of the noise parameter = 3. The modified velocity-Verlet algorithm is used for the integration over time with the time step of 0.04 [34]. For spherical particles packed into a body-centered-cubic lattice and for lamellae, the box size 60x60x60 is taken. For hexagonally packed cylinders, the box size is 60x x60. The polymer chains consist of 1000 beads each, the mean square distance between neighboring beads in a chain after equilibration is approximately equal to 1.03. Since the values of a and cut-off distance r c are close to each other, below we use in the notation for the computer simulations the value a as a unit length as if a = r c . All considered simulation boxes contain an integer number of structure periods to satisfy the periodic boundary conditions. After the equilibration, the types A and B are assigned to monomer units in accordance with their spatial positions. The monomer units in the yellow domains (lamellae, cylinders, and spheres in figure 1) are of type A and the others in the blue domains are of type B. Then, the numbers of blocks are found directly for all possible block lengths. The number and mass average block lengths are calculated and the block length distributions are plotted in comparison with the theoretical predictions.
3. Results and discussion
The blocks of type A and blocks of type B form alternating flat layers of thickness L A and L B , respectively (figure 1a). Statistical characteristics of the block distribution of type, =A,B, are determined by the layer thickness L and the monomer unit size a . The solution of the equation (1) in a layer 0 x L can be represented as a product ),|,()|,()|,()|,( ykyGzkzGxkxGkG rr where the unidimensional Green function )|,( xkxG obeys the equation xGakG . (17) Its solution is sinsin6exp2|, j L xjLxjkL ajLxkxG (18) Using equations (6) and (18), the block length distribution can be found: ,...5,3,1 02222 ~~sin~6exp~32 j LxjkLjjLkP , (19) where aLL ~ , axx ~ . In computer simulations, all macromolecules were scanned from one end to another and the block lengths and their numbers n ( k ) were found in the system after the equilibration. The total number of blocks of type is equal to k knn )( (tot) and the block length distribution can be calculated as (tot) )()(~ nknkP . If k max is the length of the longest block in the system, then kP at k > k max . If k exceeds several tens the dependence )(~ kP fluctuates remarkably, and for large k some block lengths are not found in the system. To find a smooth dependence at large values of k , we calculate the averaged fraction of blocks in an interval of length k as )(~1)(~ kk ki iPkkP , (20) where kkk . The block length distribution P ( k ) (Eq. (19)) obtained via solving the diffusion-type equation is continuous and it changes sharply for quite small k ( k < 3). For comparing this dependence with the discrete one obtained in computer simulations, we transform the expression (19) into a discrete form )(~ kP as )()(~ kPkdkP kk , (21) for which the normalization condition is fulfilled for the sum: )(~ k kP kPdk . (22) The block length distribution in a discrete form is ,...3,1 ~6)1(~60 ~~sin14)(~ j kLjkLj eeLxjjkP , k = 1,2,… (23) The probability distributions kP ~ hereafter are calculated numerically using the original program code in Free Pascal. The difference in the values kP ~ and P ( k ) given by the expressions (19) and (21) is distinguishable only at k < 5. The values of P ( k ) are calculated for the block length from k = 1 to )( max theor k corresponding to the probability P ( k ) values 10 . The accuracy of calculations is controlled by the value of the sum in the normalization condition (22), which is demanded to deviate from unity not larger than in 1-2%. The required number of terms in the sum (23) is less than ten at k > 10, whereas one hundred or more terms are necessary at k = 1. In figure 2 the block length distributions for lamellae of thickness L = L A = 5 a and 15 a are presented in the double logarithmic coordinates. The simulated values )(~ A kP are shown by markers at k < 30 and the averaged values )(~ A(av) kP are shown for larger block lengths. We increase the gap k from k = 3 at k = 30-50 to several tens for the largest observed blocks. The slope of the dashed line describing the asymptotic dependence ~ k is very close to the slope of the calculated curves and of the marker sequences at quite small k . The maximum deviation from these initial linear dependences is observed at the block lengths k aLk shown by the dotted lines, for which the second absorbing plane becomes reachable for a chain trajectory. At k > k *, the numbers of trajectories in the shape of loops and bridges become approximately equal o each other, where “loops” and “bridges” are trajectories absorbed at the planes x = 0 and x = L , respectively [15]. Figure 2. Block length distributions )(~ A kP for blocks in the lamellae of thickness L = L A = 5 a (blue curve and markers) and 15 a (black curve and markers). The curves are calculated from Eq. (23) at x = 0.5 a . The markers describe the results of the computer simulations. The line describing the dependence ~ k is plotted to guide the eye. The vertical dotted lines correspond to the values of aLk . In figure 3, the number average block size A nn kk and the dispersity Ð = Ð A that characterizes the relative width of the distribution )(~ A kP are presented dependently on the lamellar thickness L = L A . These values are calculated for the discrete distributions )(~ A kP as ~ kn kPkk , ~1 knw kPkkk , Ð nw kk (24) The dependences of n k and Ð on the lamella thickness are practically linear in the numerical calculations. For the simulation data, the deviation from the linear law is observed at L / a = 15, where the presence of chain ends probably influences the distribution width. Figure 3. Number average block length A nn kk (thin solid curve) and dispersity Ð = Ð A (thick solid curve and markers) for blocks in the lamella vs lamellar thickness L = L A . Solid curves and markers describe the results of the numerical calculations and the computer simulations, respectively. Blocks of type A are located in the hexagonally arranged cylinders of radius R surrounded by the continuous B domain (figure 1b). The distance between the axes of nearest cylinders is equal to d . The analytical solution of Eq. (1) for the Green function G A ( r , k | r ) in a cylinder is well known [31], however, the solution for an infinite domain with excluded cylinders has not been reported yet. Any block begins near a cylinder surface that crosses a polymer chain. Any A block is localized within one cylinder. For a B block with the beginning near a certain cylinder, we approximate the presence of other cylinders by a coaxial cylindrical surface of radius r out (figure 4). Figure 4. The cross-section of the melt with hexagonally arranged cylinders of radius R with A blocks surrounded by the external domain with B blocks. The distance between the axes of neighboring cylinders is equal to d (the structure period). For the B blocks beginning near the central cylinder, the presence of the other cylinders is described in the theoretical model by the cylindrical surface of radius r out . The three-dimensional Green function of sequences can be represented as a product )|,(,|,,)|,( zkzGrkrGkG rr , where the Green function ,|,, rkrG obeys the equation GrrGrrrakG (25) The function G A2 is defined for a circle at 0 r R and 0 , the boundary condition is rkRG . The function G B2 is defined for a ring at R r r out and 0 , the boundary condition is rkrGrkRG . The initial condition for the both functions is )()(1 rrrG k . Since the block length distribution does not depend on the choice of the orientations of the coordinate axes, the convenient value = 0 can be used. The Green function G A2 for a circle has the form j n njnjnnjnnjnn kRaRrJRrJJ nbRG (26) where J n is the Bessel function of the n -th order, b n = 1 for n = 0 and b n = 2 for n )( nj is the j -th positive solution of the equation J n ( ) = 0. The values of these solutions for the Bessel function of the zeroth order are )0(1 = 2.40, )0(2 = 5.52, )0(3 = 8.65, )0(4 = 11.79, )0( 1 j )0( j . r out R d
A B The block length distribution P A ( k ) for blocks in cylinders can be found using Eqs. (5) and (24) as A20 20A )( GddrrkkP R . Since the terms with n > 0 are equal to zero, the distribution P A ( k ) can be written as kRfRrJJRkP jjj jjj ~6 )(exp)(~~)(~3 1)( , (27) where
10 00 xxJdxf , aRR ~ , arr ~ . Using the relations between the Bessel functions of the 0th and the 1 st order, P A ( k ) can be written as kRRrJJRkP jj jjj ~6 )(exp~~~3 1)( (28) Then, the probability distribution in a discrete form is kRkRJ RrJkP jjj jj j ~6 )(exp)1(~6 )(exp~~2)(~ . (29) The Bessel functions were calculated numerically using the necessary numbers of terms in the sums n nn xnxJ , n nn xnnxxJ for x < 37. For larger arguments, the asymptotic form )4sin(81)4cos(2 xxxxxJ was taken, xJxJ . The Green function G B2 for B blocks in a ring has the form j n njnjnnjnnjnnjn njnnjn kRRrfRrfcJJ ncJbRG , (30) where cNJcJNcf nnnnn , , N n is the Bessel function of the second kind (the Neumann function) of the n -th order, )( ~ nj is the j -th positive solution of the equation out cf n , Rrc outout . Using equations (5) and (30), the probability distribution P B ( k ) is calculated as B220B )( GddrrkkP out rR : ~6 )~(exp),~(~~,~~~ ~)~(~12)( j joutjjoutjj outjjB kRcfRrfcJJ cJRkP , (31) where ),(),( xfxdxcf c . Here, t he Neumann function of the zeroth order can be calculated as [35] ),()(2log2)( h00 xfxJxxN mmm m hxmxf
21 2h , where = 0.5772157… denotes the Euler constant, h m = 1 +2 +…+ m . The function ),( B cf can be rewritten as ),,(1 )(21ln2),( h10 h110010B cfJ fJccJJcJcccJJcf (32) where )1(2out21 2h1 mmm mm cxmm hcf . The values of ),( B cf were calculated using the expression (32) at c < 37, at larger arguments the function ),( cf was integrated numerically. The probability distribution P B ( k ) in a discrete form is kRkR cfRrfcJJ cJkP jjj outjjoutjj outjjB ~6 )~(exp)1(~6 )~(exp ),~(~~,~~~ ~)~(2)(~ (33) The results for the probability distributions )(~ A kP for blocks in the cylinder of radius R = 4 a and )(~ B kP for blocks in the external domain are presented in figure 5. The external domain in the computer simulations encompasses the cylinders of radius R = 4 a packed into the hexagonal structure with a period d = 15 a . In the numerical calculations, the probability distribution )(~ B kP is presented for two cylindrical layers between R = 4 a and two outer radii r out = 15 a (thin blue curve) and r out = 13.1 a (thick blue curve). The second curve provides the best fit for the slope of the dependence of )(~ B kP on k in the inset graph. Figure 5. Block length distributions )(~ A kP for A blocks in the cylinder of radius R = 4 a and )(~ B kP for B blocks in the external domain. Markers represent the results of the computer simulations for the cylinder (triangles) and for the external domain (rhombuses) of the hexagonal structure with the period d = 15 a . Solid curves describe the results of the numerical calculations for the cylinder (thick black curve, R r = 0.5 a ) and for the cylindrical layer with R / a = 4, r R = 0.5 a , r out / a = 15 (thin blue curve) and 13.1 (thick blue curve). The inset represents the same dependences in the plot with the logarithmic scale only at the vertical axis. In figure 6, the number average block length A n k and the dispersity Ð A for A blocks in the cylinders are plotted dependently on the cylinder radius R . The dependences are practically linear both in the theory and computer simulations. In figure 7, the number average block length B n k and the dispersity Ð B for B blocks in the external domain are plotted dependently on the structure period d in the computer simulations or on the outer radius r out in the numerical calculations at the fixed cylinder radius R . The same parameters are presented in figure 8 depending on the cylinder radius R at the fixed sizes of the external domain. Figure 6. Number average block length A n k and dispersity Ð of blocks A in the cylinders vs cylinder radius R / a . Solid curves with markers describe the results of the numerical calculations for A n k (thin curve with triangles) and Ð (thick curve with circles). The separate markers describe the results of the computer simulations for A n k (blue triangles) and Ð (red squares). Figure 7. Dependence of the number average block length B n k (a) and dispersity Ð B (b) of B blocks in the external domain on spatial period d / a (thin curves with markers) at the values of cylinder radius R / a = 3 (1), 4 (2), 5 (3), and 6 (4) in the computer simulations and on outer radius r out (thick black curve) at the radius R = 5 a in the numerical calculations. a b Figure 8. Dependence of the number average block length B n k (a) and dispersity Ð B (b) of B blocks in the external domain on cylinder radius R / a in the computer simulations (thin curves with markers) at the values of spatial period d / a = 10 (1), 12 (2), 15 (3), 20 (4), 30 (5) and in the numerical calculations (thick black curves) at the outer radius r out =12 a . Balls of radius R s are arranged in a body-centered cubic lattice, those balls consist of A blocks and the surrounding domain consists of B blocks (figure 1c). The distance between the centers of nearest balls is equal to d s . The solution of Eq. (1) for sequences of A monomer units in a ball can be written in a usual way [31], whereas the form of the solution for an infinite region with excluded balls is not known. For a B block beginning near a certain ball, the presence of other balls is approximated by an absorbing spherical surface of radius r out with the same center. The Green function of -type sequences obeys Eq. (1) that can be written in the spherical coordinates for the function G ( r , , , k | r , , ) as sin1sinsin16 GGrGrrrakG (34) The function G A is defined for a ball at 0 r R s , 0 , and 0 ; the boundary condition is sA Rr G . The function G B is defined for a spherical layer at R s r r out , 0 , and 0 ; the boundary condition is BsB out rrRr GG . The initial condition for the both functions is )()()(1 rrrG k . The convenient orientations of the coordinate axes can be chosen to provide the initial values = 0 , = 0. The Green function G A for the ball has the form a b kR RrJRrJJ PnrrRG njj n njnnjnnjn n
2s )21(1 0 s)21(21s0)21(212)21(2102sA ~61exp)(cos)12(2 1 (35) where P n ( x ) is the Legendre polynomial, )21( nj is the j -th positive solution of the equation n J , aRR ss ~ . Since sin2 J , these solutions at n = 0 are j j )21( . The block length distribution P A ( k ) for blocks in the ball can be found as A0 2002A s sin)( GdddrrkkP R , where only the terms with n = 0 give non-zero contributions: ~6expsin)1(~~3 j j kRjRrjjrRkP (36) In a discrete form, ~6exp)1(~6expsin)1(~~2~ j j kRjkRjRrjjrRkP (37) The solution of Eq. (34) for the function G B in a spherical layer at = 0, = 0 can be written in the form kRRrfRrf RrJJ RrJPnrrRG njsnjnsnjnj n snjnnjn snjnnjn
2s 2)21()21(210)21(211 0 out)21(2 21)21(2 21 out)21(2 212)21(02sB ~ )~(61exp),~(),~( )~()~( )~()~)((cos)12(8 (38) where cNJcJNcf nnnnn , , )21( ~ nj is the j -th positive solution of the equation out21 cf n , s Rrc outout . The block length distribution P B ( k ) can be calculated as B2002B out sin)(
GdddrrkkP rR s . The integrals of P n (cos ) over are not equal to zero only for n = 0, therefore, only the terms with P (cos ) = 1 should be taken. The equation out21 cf is equivalent to out c . Introducing the notations R = r out R s , r = r R s , aRR ~ , the block length distribution can be written in the form kRjRrRrjjRr aRkP j sjs
221 out020 2B ~6 )(exp)1(1sin3)( , (39) where the value
RjR sj )21( ~ is used. In a discrete form, kRjkRjRrRrjjrRkP sjjs ~)(61exp)1(~6 )(exp)1(1sin12)( (40) The block length distributions )(~ A kP for A blocks in the cylinder and in the ball of the same radius 5 a are presented in figure 9. The distributions found in the computer simulations are very close to each other at k < ( R / a ) , whereas shorter blocks prevail in the ball at larger k in comparison with the cylinder. The block length distributions )(~ B kP of B blocks in the external domains of different sizes are plotted in figure 10 at R s = 5 a . The solid curves are plotted for the numerical data at the values of r out that provide the best fit for the slope of the dependences obtained in the computer simulations at large k . Figure 9. Block length distributions )(~ A kP of A blocks in the cylinder of radius R = 5 a (green triangles) and in the ball of radius R s = 5 a (blue rhombuses and solid curve). The markers represent the results of the computer simulations, the solid curve describes the numerical calculations, R r = 0.5 a . Figure 10. Block length distribution )(~ B kP of B blocks in the external domain surrounding a ball of radius R s = 5 a at d s / a = 12 (1), 15 (2), and 20 (3) in the computer simulations (markers) and at r out / a = 8 (1), 11.5 (2), and 18 (3) in the numerical calculations (black solid curves, r R = 0.5 a ). Figure 11. Number average block length A n k and dispersity Ð A of blocks A in a ball vs ball radius R s / a . Solid curves describe the results of the numerical calculations for A n k (thin curve) and Ð (thick curve). Markers describe the results of the computer simulations for A n k (blue rhombuses) and Ð (red squares). In figure 11, the number average block length A n k and the dispersity Ð A for A blocks in a ball are plotted dependently on its radius R s , and they grow approximately linearly with the radius as well as for a cylinder. In figures 12 and 13, the dependences of the number average block length B n k and the dispersity Ð B for B blocks in the external domain are plotted vs the structure period d (or the outer radius r out ) and the ball radius R . Figure 12. Dependence of the number average block length B n k (a) and dispersity Ð B (b) of B blocks in the external domain on spatial period d s / a (thin curves with markers) at the values of ball radius R s / a = 3 (1), 5 (2), and 6 (3) in the computer simulations and on the outer radius r out (thick black curves) at the radius R s = 3 a in the numerical calculations. Figure 13. Dependence of the number average block length B n k (a) and dispersity Ð B (b) of B blocks in the external domain on ball radius R s / a in the computer simulations (thin curves with markers) at the values of spatial period d s / a = 12 (1), 15 (2), 20 (3) and on the outer radius r out / a in the numerical calculations (thick black curves) at the outer radius r out =10 a . a b a b In the previous three subsections the formulae for the block length distributions are derived for three considered structure types (figure 1), and the dependences of the probabilities )(~ kP on the block length k and of the distribution parameters n k , Ð on the domain sizes are plotted. All the probabilities )(~ kP (figures 2, 5, 9, and 10) obviously obey the asymptotic law ~ k at the moderate values of block length (the asymptotic line is shown in figure 2). At very large k , the probabilities decrease exponentially in the theory and practically by the same way in the computer simulations. The block length distributions )(~ A kP obtained with the computer simulations and with the numerical calculations at the same values of the domain size (lamella thickness L A , cylinder radius R , or ball radius R s ) are close to each other (figures 2, 5, 9). The distribution for a cylinder is wider than for a ball of the same radius 5 a (figure 9), since A blocks are bounded in all three dimensions in the ball, whereas there are no boundaries along the cylinder axes. The ratio of the number average values of block lengths for a cylinder and a ball of radius 5 a is approximately equal to A(ball)A(cyl) nn kk both in the numerical calculations and in the computer simulations (figures 6 and 11) despite the values of the lengths A n k themselves are different. The probability distributions )(~ A kP for lamellae (figure 2) and )(~ B kP for the external regions of cylinders and balls (figures 5 and 10) can demonstrate non-monotonic change of the slope in the double logarithmic scale. The hump at the intermediate values of block length k corresponds to the critical value k *, for which a block with k > k * can reach an additional absorbing boundary. That value can be estimated as aLk for lamellae or aR for cylindrical and spherical layers of thickness R in the numerical calculations. In the computer simulations, such hump is less pronounced and can be related to k * = aR , where min R is the smallest distance between the boundaries of neighbor cylinders or spheres. At large k , the theoretical probability distributions )(~ A kP and )(~ B kP for the different structure types are described by the asymptotes, which are the first terms ( j = 1) in the expansions (19), (27), (30), (34), and (37). The slope of the linear dependence of )(~ln kP on k (the inset in figure 5) is inversely proportional to the square of the characteristic spatial size of the domain. The external domains of cylinders and balls (figure 1b,c) are approximated in the theoretical model by the symmetrical domains in the shape of cylindrical and spherical layers, and the width of those layers can be taken as the characteristic domain size. In the computer simulations, the probability distributions )(~ B kP for the external domains are found directly. To estimate the effective outer radius )eff(out r for such domains, the slope of the dependences of )(~ln kP on k is fitted by the corresponding theoretical solutions. For the structure of the cylinders of radius R = 4 a with the period d = 15 a , the best fit is provided by the value r out = 13.1 a (figure 5), which lies between d R and d . For the structure of balls, three dependences are fitted for the different values of d (figure 10), all of them are within the range from d s R s to d s , and we see that )eff(out r is not directly proportional to d . In the limit of large values of period ( d / R , d s / R s >>1), the effective value of )eff(out r can be approximated with the power-law dependence. Consider a ball in the center of the first coordination sphere of radius d s . The fraction s of the sphere domain covered by other balls is proportional to ~ dR , the number of the coordination spheres in the ball of radius )eff(out r is equal to s)eff(out drm . Assuming that the fraction s is approximately the same for all coordination spheres, the solid angle covered by other balls from the point of view in the center of a chosen ball is equal to m s . Let )eff(out r be a distance at which any random trajectory will be finally absorbed, then m s
1. For the system of cylinders, the first coordination cylinder of radius d should be considered, with the surface fraction ~ dR covered by other cylinders. Similar reasons lead to the condition m
1, where drm )eff(out . Then, the effective outer radius in the limit of large d Rdr ~ (cylinders), ~ Rdr (balls). (39) The dependences of number average block lengths n k and dispersities Ð on the domain sizes are also presented for all the structure types (figures 3, 6-8, 11-13). The solid curves describe the results of the analytical calculations and the markers describe the results of the computer simulations. For a lamella (figure 3), a cylinder (figure 6), and a ball (figure 11), the average block length A n k and the dispersity Ð grow practically linearly with the size of the corresponding domain ( L , R , and R s , respectively) both in the theory and in the computer simulations. The linear law (12) and (13) follows from the Levy-flight statistics P ( k ) ~ 1/ k at the moderate values of k (9) together with the exponential decrease (10) beginning from the unique characteristic size k* ( d A / a ) . However, the values of A n k found in the computer simulations are remarkably larger than those calculated numerically (see, for example, figure 6 for cylinders). This difference can be explained by the peculiarities of the analytical solution for the probability distribution P A ( k ) at small values of block length k . Chain trajectories are described by the Green function (Eq. (1)) as continuous ones, whereas a set of monomer unit coordinates is discrete. It can be expected that the continuous description is not adequate at small k . Indeed, at k = 1 the probability A P in theory and 0.32 in the computer simulations for the ball of radius R = 4 a . Since the fraction of short blocks of size k = 1 and 2 is overestimated in the analytical solution, the fraction of other blocks is underestimated thus leading to smaller values of the number average block size A n k in the analytical solution than in the computer simulations. For the values of the mass average block size w k and dispersity Ð , the largest contribution is made by the longest blocks. In the computer simulations, a block trajectory can be interrupted by a chain end before it crosses the interface. A lack of long blocks in the computer simulations in comparison with the analytical solution leads to somewhat smaller values of the dispersity (figures 3, 6, and 11). For B blocks, the dependences of number average block length B n k on the outer radius r out and inner radius ( R or R s ) obtained in the analytical calculations are much smoother than the dependences of B n k on the size of the external domain (spatial period or radius) obtained in the computer simulations (figures 12a, 13a, 15a, and 16a). In the analytical calculations, the dependences are practically linear in the presented cases. In the computer simulations, B n k grows faster than linearly with the structure period d (or d s ) under the fixed cylinder (or ball) radius. For a block beginning near a given cylinder, the nearest points of the surfaces of the other cylinders move away with increasing d , and simultaneously other cylinders become rarefied in space. In the result, the effective outer radius should change faster than linearly with d . And the same conclusion is valid for the structure of balls arranged in the lattice. The dispersity Ð B changes practically linear with the period d ( d s ) and with the radius R ( R s ) at moderate values of the thickness d R in the computer simulations. For larger values, the dispersity sharply decreases mainly because of a decrease in the values of mass average block length B w k , which are sensitive to the lack of long blocks. For example, there is a certain fraction of the homopolymer B chains of length 1000 in the simulation box for the balls of radius R s = 3 a with the structure period d s = 20 a . Therefore, in the limit N a noticeable fraction of longer blocks is expected to be present in this system. The value of Ð B calculated from the expressions (24) and (33) changes approximately linearly with r out and R (thick curves in figures 7b and 8b) at the comparable values of the cylinder radius R and the thickness r out R . For the length distribution (40) of blocks in the spherical layer with the inner radius R s and outer radius r out , the dependence is also close to linear (figure 12b and 13b). The dispersity is quite large in all considered cases and, besides, it is predicted to grow with the domain size (estimation (13)). That is, the dispersities can considerably exceed those for the polymers synthesized via radical polymerization, which is characterized by the most probable (Flory) distribution with the dispersity equal to 2. Despite a very broad distribution of block lengths, the polymers obtained via the conformation-dependent sequence design are able to form stable segregated structures. The stability of globules with a hydrophobic core and a polar shell formed by protein-like copolymers was demonstrated in the computer simulations and in the theory [7, 8] and also experimentally [9, 10]. The formation of lamellar structure in the melt of multiblock copolymers prepared via the design with a lamellar pattern was also observed in the computer simulations [15]. After switching on the repulsive interactions between monomer units of different types, the lamellae with the period considerably larger than that of the initial pattern were formed. It can be expected that other patterns can also be reproduced for the other types of the designed block copolymers in melt. It is worth to note that very short blocks can be located in the domains of “alien” type or at the interface and, therefore, they hinder the microphase separation in a melt. The effect of short blocks on the microphase separation was studied theoretically for the melt of the random multiblock copolymer of special type [36]. Thus, for better reproduction of a pattern used in the conformation-dependent design, it is desirable to prevent the formation of very short blocks during the modification.
4. Conclusion
In the present paper, the statistical characteristics of random multiblock copolymers of the special type have been studied. The type A or B is assigned to the monomer units in accordance with their spatial positions in the homopolymer melt, where the periodic pattern is introduced. The monomer unit sequences designed in such a way depend on the pattern type and the domain sizes (figure 1). The designed random multiblock copolymers are able to reproduce, for example, the “parent” lamellar structure in the melt with repulsive interactions of A and B monomer units, as was demonstrated previously in the computer simulations [15]. We have been considered together the several structure types (alternating lamellae, hexagonally arranged cylinders, and balls packed in the body-centered cubic lattice), which are typical of microphase separation in block copolymer melts. Such bulk modification of a polymer melt still has not been performed experimentally. As a way to realize the process, we can suggest to initiate a photochemical reaction in a sample creating a light 3D-pattern by optical methods. For the modified macromolecules, the block length distributions and the average block lengths have been found both in the computer simulations of the melt equilibrated via the DPD method and theoretically using the description of random walk statistics of polymer chains with the diffusion-type equations. When moving along a polymer chain, the monomer unit type changes at the interfaces. Simultaneously, the polymer block is terminated, as described by the absorbing boundaries or Dirichlet-type boundary conditions. For the external domain of balls and of cylinders, the approximation of a spherical layer and a cylindrical one, respectively, have been used in the theory (figure 4). At moderate length k of blocks, the Levy flight-type distribution of block length P ( k ) ~ 1/ k is observed for all the pattern types. This dependence describes the asymptote of the block length distribution for the blocks starting near an absorbing plane in a semi-space. Therefore, this asymptote should be valid for any pattern in 3D until the initial surface curvature reveals itself and the other boundaries become attainable, that is, before the characteristic spatial size ak of a block reaches the smallest characteristic size of the domain. The general shape of the block length distributions obtained in the computer simulations and in the theoretical calculations is very similar for the blocks in lamellae, cylinders, and balls. For the external domain of balls and of cylinders, those distributions are not so similar. However, the slope of the dependence found in the computer simulations can be well fitted at large block lengths by a proper choice of the thickness of the spherical or cylindrical layer in the theoretical calculations. The asymptote exp~)( dkakP , = const, can be used at large k , where the asymptotic length d as nonlinearly depends on two parameters, R and d (or R s and d s ), for the external domains of cylinders and balls. It can be expected that for more complex lattices such asymptotic length can depend on several geometric parameters that determine the size and arrangement of the domains. For the symmetric domains of cylindrical and spherical layer, the layer thickness is the only relevant parameter. The average block length and the dispersity depend linearly on the characteristic domain size if this size is unique. In a multi-scaled domain, such dependence can be nonlinear. The numerical values of the average block lengths calculated from the analytical solutions of the diffusion-type equations are somewhat underestimated in comparison with the results of the computer simulations because of the excess number of very short blocks, for which the continuous theoretical description is not very adequate. On the other side, the block length distribution in finite polymer chains in the computer simulations does not describe the statistics of infinite random walks for patterns with a large spatial period. In that case, the effect of chain ends, which terminate the blocks, becomes important for the chains of spatial size comparable with the structure period. Therefore, both methods complement each other and can provide rich information about the statistics of random walks in a patterned media. Acknowledgment
The research was financially supported by the Russian Foundation for Basic Research (project no. 19-03-00988).
References [1] Runnels C M, Lanier K A, Williams J K, Bowman J C, Petrov A S, Hud N V, and Williams L D 2018 Folding, Assembly, and Persistence: The Essential Nature and Origins of Biopolymers.
J. Mol. Evol. Bioinformatics Introduction to Protein Structure Prediction: Methods and Algorithms (ed. by Rangwala H and Karypis G)
Reviews of Modern Physics Phys. Rev. Lett. , 3456 [6] Govorun E N, Ivanov V A, Khokhlov A R, Khalatur P G, Borovinsky A L, Grosberg A Yu 2001 Primary sequences of proteinlike copolymers: Levy-flight-type long-range correlations. Physical Review E R40903 [7] Govorun E N, Khokhlov A R, and Semenov A N 2003 Stability of dense hydrophobic-polar copolymer globules: Regular, random and designed sequences.
Eur. Phys. J. E Phys. Rev. E Adv. Polym. Sci. and studies of N-vinylcaprolactam/N-vinylimidazole copolymers that exhibit the “proteinlike” behavior in aqueous media. Macromolecules Eur. Phys. J. E J. Chem. Phys.
Phys. Rev. E Macromolecules
42 2843-2853 [15] Govorun E N, Gavrilov A A, Chertovich A V 2015 Multiblock copolymers prepared by patterned modification: Analytical theory and computer simulations.
J. Chem. Phys.
Biomacromolecules Eur. Polym. J.
Beilstein J. Org. Chem. ChemRxiv doi.org/10.26434/chemrxiv.12895754.v1 [20] Gody G, Maschmeyer T, Zetterlund P B, and Perrier S 2013 Rapid and quantitative one-pot synthesis of sequence-controlled polymers by radical polymerization.
Nature Communications Statistical Physics of Macromolecules (Woodbury, NY: AIP) [22] Feller W 1971
An Introduction to Probability Theory and Its Applications , vol. II (3rd ed.) (New York: John Wiley & Sons) [23] Redner S 2001 A Guide to First-Passage Processes (Cambridge University Press) [24] Bray A J, Majumdar S N, and Schehr G 2013 Persistence and First-Passage Properties in Non-equilibrium Systems.
Advances in Physics
62 225-361 [25] Ben-Naim E and Krapivsky P L 2010 First-passage exponents of multiple random walks. J. Phys. A: Math. Theor. Risks
117 [27] Simon D and Derrida B 2008 Quasi-Stationary Regime of a Branching Random Walk in Presence of an Absorbing Wall.
Journal of Statistical Physics
J Korean Phys. Soc.
Numer. Math. The Fractal Geometry of Nature (Updated and augm. ed.) (New York: W. H. Freeman) [31] Budak B M, Samarskii A A, and A N Tikhonov A N 1964
A collection of problems on mathematical physics (Pergamon Press) [32] Hoogerbrugge P J and Koelman J M V A 1992
Europhys. Lett. J. Chem. Phys.
Computer Simulation of Liquids (Oxford: Clarendon Press). [35] Watson G N 1944
A treatise on the theory of Bessel functions (Cambridge University Press) [36] Govorun E N and Chertovich A V 2017 Microphase Separation in Random Multiblock Copolymers
J. Chem. Phys.146