Charged designed copolymers in the presence of multivalent counterions: a molecular dynamics study

We present the results of molecular dynamics simulations of charged proteinlike hydrophobic–hydrophilic (HP) copolymers in a dilute salt-free solution with multivalent counterions under poor solvent conditions. The primary sequence of these copolymers is constructed such that they can self-assemble into a segregated core–shell microstructure thus resembling some of the basic properties of globular proteins. The results are compared with those obtained earlier for the system containing monovalent counterions. The processes of coil-to-globule transition, aggregation of polyions, and counterion condensation are studied in detail as a function of temperature. The main attention is paid to the influence of the counterions of different charge on the aggregation of the copolymers. It is found that multivalent counterions considerably increase the aggregation of the chains, which form in their presence the finite-size aggregates built up from several polyions. However, the striking feature of the aggregation is that this process does not appear to lead to macroscopic phase separation. The mechanism that prevents the phase separation is discussed.


Introduction
It is well known that many of the solution properties of biopolymers (DNA, proteins) are due to a complex interplay between short-range hydrophobic effects and long-range Coulomb interactions. In particular, the balance of these interactions determines the solubility and equilibrium conformation of proteins in water [1]. To a large extent, these properties are also dependent on the unique primary structure of the protein chain [1]. This fact motivates the experimental and theoretical studies directed to the design of synthetic polyelectrolytes with specific primary sequences resembling those of biomacromolecules.
In [2]- [4], an approach to the conformation-dependent sequence design was realized for the so-called 'proteinlike' HP copolymers consisting of hydrophobic (H) and hydrophilic (P) monomer units. The primary sequence of these copolymers is constructed such that they can self-assemble into a segregated core-shell microstructure thus resembling some of the basic properties of globular proteins, namely, their solubility in water. Both theoretical and computer simulation studies showed that the behaviour of proteinlike copolymers is essentially different from that of the random and random-block copolymers with the same HP composition and the same degree of blockiness [2]- [6]. In particular, Govorun et al [6] found that the proteinlike primary sequence demonstrates the specific long-range correlations, which can be described by the statistics of the Lévy-flight type [7]. All of these studies focused on uncharged proteinlike copolymers whereas many of the biopolymers are polyelectrolytes.
Recently, we have reported the results of molecular dynamics simulations of charged proteinlike HP copolymers with a fixed charge distribution under poor solvent conditions in the presence of monovalent counterions [8]. Three temperature regimes characterized by a different behaviour of the Coulomb energy, chain sizes and pair correlation functions were found. In the high-temperature regime, the polyions had extended coil-like conformations with many hydrophobic blobs while counterions were distributed almost uniformly in solution. As the temperature was decreased, the chains formed globules with electroneutral parts located in the globular core and charged sections forming the envelope of the core and buffering it from the polar solvent. In this intermediate-temperature regime, we observed a stable solution of non-aggregating polymer globules, which are well separated from each other and form an array of colloid-like particles. A further decrease in temperature led to the condensation of most of the counterions and to the aggregation of the individual globules which, nevertheless, maintained their morphological integrity in finite-size aggregates. A solution of globules has also been studied in the molecular dynamics simulations of Kremer and coworkers [9,10], who considered highly charged regular copolymers in a poor solvent in the presence of monovalent counterions. Recently, Zherenkova et al [11] have performed rather extensive studies of the static structure and self-assembling of charged quasi-random copolymer solutions on the basis of the integral equation (PRISM) theory. To gain insight into the influence of the primary structure of copolymers on their structural behaviour, two types of copolymers were considered: proteinlike and random-block copolymers. It was shown that charged proteinlike HP copolymers are more prone to self-organization in a poor solvent than their random-block counterparts. The results obtained for solutions containing polyions and monovalent counterions have also shown that the differences between two copolymers with different primary structure in the cases where the electrostatic interactions play the main role are much smaller than the effects driven by short-range attractive interactions between hydrophobic chain segments [11]. It may well be, however, that the situation will be different when multivalent counterions are present in the solution.
It is known that the interaction of polyelectrolytes with multivalent counterions can be a critical factor both in many industrial applications and in biological functions. The presence of multivalent counterions in solution can have a dramatic effect on the aggregation of the polyions. In particular, strongly charged polyelectrolyte chains can condense at a critical concentration of multivalent ions. Usually, this process is followed by a macroscopic phase separation [12]. But in recent years it has been shown that some biopolymers can condense in the presence of multivalent ions into soluble aggregates of a well-defined size. Such a phenomenon was observed for doublestranded long DNA dilute solutions and for short DNA fragments [13]- [18], for F-actin [19], for tobacco mosaic and fd viruses [20], etc. The existence of finite-size thermodynamically stable multichain aggregates in the dilute solutions of associating polyelectrolytes has also been predicted theoretically in a series of papers [21]- [24].
In this paper we study the effect of multivalent ions on the processes occurring in dilute salt-free solutions of charged proteinlike copolymers. Our main aim is to study conformations of such copolymers, as well as their stability towards large-scale aggregation under poor solvent conditions.

Model and simulation technique
The simulated system is a dilute solution of N -unit polyelectrolyte chains with a fixed charge distribution and oppositely charged small mobile counterions. There is no additional salt in the system. The solvent is replaced by a dielectric background and a heat bath. We consider a simple model of a flexible copolymer chain that involves only two types of monomer units, H and P. The P units are assumed to be negatively charged, while the H units are electroneutral. The HP composition is fixed at 1 : 1. We employ a continuum space (bead-rod) model of a polymer chain. The time evolution of the system is determined by Newton's equations that are solved using the method of molecular dynamics (MD). Monomer units (chain beads) are linked by rigid bonds of a fixed length to form a linear HP copolymer chain of length N .
The two-particle potential energy of the system is given by where r ij = |r i − r j | is the distance between interacting particles i and j, N is the total number of particles (monomer units plus counterions) in the system, u ev takes into account excluded volume, u a characterizes attractive interactions between monomer units, u C is the energy of electrostatic interaction between charged particles (P monomer units and/or counterions). Excluded volume between any non-bonded particles (i.e., not connected by a bond) is included via a repulsive Lennard-Jones (RLJ) potential where σ = ε = 1 and r 0 = 2 1/6 σ is a cutoff distance. The parameter ε entering this equation controls the energy scale, whereas σ determines the interaction length between particles. In addition to the excluded volume potential (2), the intramolecular and intermolecular interactions between non-bonded chain beads are modelled by a Yukawa-type potential, for which we use the following form [8]: In the model, this potential describes short-range attractive intra-and interchain interactions. The parameter ε αβ (= ε HH , ε PP , ε HP ) sets the depth of the minimum of the attraction and r c = 2.8σ is the cutoff distance for attractive interactions. The characteristic energy of H-H, H-P, and P-P interactions is fixed at ε HH = ε HP = ε PP = 2. For ε αβ = 0, there is no attraction between chain sites. For a system of N particles with charges q i at positions r i in a cubic simulation cell of length L box with periodic boundary conditions, the long-range Coulomb energy is given by where r ij = r i − r j , q i = ez i is the charge of a given particle, e the elementary charge, ε 0 denotes the permittivity of free space, ε r the relative permittivity of the medium, and the sum over on polyions and counterions are equal to −1and +z, respectively. We perform calculations at z = 1, 2 and 4. In each case, we consider the polyelectrolyte solution containing only one type of the counterions, monovalent, divalent and tetravalent ones. Since in our simulation all charges are treated explicitly, electrostatic interactions are described by the electrostatic potentials (4) with a dielectric permittivity ε r = 1. Additional computational details and parameters used in the simulation are given in the appendix. Explicitly, no solvent molecules are included in the simulation; i.e. the solvent is represented by a dielectric continuum. In order to simulate solvation effects and the time evolution of the solution in contact with a heat bath of temperature T , we augment the equations of motion by a frictional term and a Langevin uncorrelated noise term [25,26]. In the MD simulations, the equations of motion of a polymer system in the presence of fixed constraints (the bead-rod model) are with the constraints of fixed bond lengths where n is the number of chains in the system, m i = 1 the mass of the chain bead or counterion i, −∇ i U(r) represents the total non-bonded forces on the particle i, b = σ the bond length, λ k s are the Lagrange multipliers, R describes the random force of the heat bath acting on particles, and ξ is the frictional (damping) coefficient which takes into account the viscosity of the solvent. The summation in (5) is performed over all the n(N − 1) bonds of the N -unit polymer chains. The term ∇ i n(N −1) k=1 λ k G k (r 1 , . . . , r nN ) represents the forces, which are due to the bond reactions. In (5), R i and ξ i are connected through the fluctuation-dissipation theorem, R αi (0) · R αj (t) = 2mξk B Tδ ij δ(t), (α = x, y, z; T is the reference temperature) which ensures that the temperature is kept constant. We take the parameters ξ i to be dependent on solventaccessible surface areas (SASA) (for more detail, see [27]). To find the values of SASA for a given configuration, we perform an analytical computation of the surface areas A i for each specified particle [28,29]. Having A i , one can define ξ i as ξ i = ξ 0 A i /A max , where A max is the maximum solvent-accessible surface area of a particle for the model under study and the reference value of ξ 0 is set to be equal to unity (in the range between overdamped regime and purely deterministic regime).
The equations of motion (5) in conjunction with equations (6) are solved iteratively using a Newtonian iteration procedure with the time step t = 0.01 σ √ m/ε. For this calculation we use the algorithm described in [30,31]. The unit of time in the simulation is defined by τ = σ √ m/ε. In the following, the reference temperature T (measured in units of ε/k B ) is considered as the main variable parameter.
In order to obtain single proteinlike HP copolymer, the following computational strategy was employed [2]- [4]: first, a self-avoiding homopolymer chain with uncharged H units is randomly generated. Then the folding of the chain is performed at ε HH = 2 and T = 1. This conformation is equilibrated for 4 × 10 5 integration time steps. As a result, we obtain a dense homopolymer globule. One half of the monomer units, having the largest SASA and simultaneously wellseparated from the centre-of-mass of the globular core, are transformed into P type; the remaining monomer units with the (minimum) values of SASA and being closer to the centre-of-mass of the globular core are considered to be of H type. Thus, the composition of the polymer sequence is constrained so that there are N /2 charged and N /2 neutral monomer units, and the primary structure of the heteropolymer is frozen. As a result, we obtain a bare proteinlike HP heteropolymer globule having a core-shell structure.
The primary structure of a two-letter HP copolymer can be characterized by its composition, by the average lengths of the H and P blocks, L H and L P , and by the specific distribution of H and P units along the chain.
The bare globule was replicated to a total number of 29. To this end, the images of the bare globule were placed on a regular simple cubic lattice spanning a periodic MD box without overlap between them and then were randomly reoriented. In subsequent calculations, each P monomer unit is assumed to be negatively charged, while H units are electroneutral. Then, positively charged counterions were distributed in the box and the system constructed in this way was relaxed for a long time.
For computational efficiency, chains with relatively short length, N = 128, were used for most of the calculations. The cubic computational box contained n = 29 polyions (total number of chain beads nN = 3712) and n c = nN /2z counterions to ensure electroneutrality of the system. The entire system of n chains plus n c counterions was enclosed in a periodic cubic box with side length L box (L box = 48; volume V = L 3 box ). The number density of the chain beads is ρ p = nN /V = 3.356 × 10 −2 and the number density of counterions is ρ c = n c /V (ρ c = 1.678 × 10 −2 at z = 1, 8.38 × 10 −3 at z = 2 and 4.19 × 10 −3 at z = 4). The value of ρ p is sufficiently low to correspond to a dilute polyelectrolyte solution [8].
The systems were simulated at different temperatures in the range from T = 0.25 to 7.0. This temperature range corresponds to a Bjerrum length l B = e 2 /4πε 0 ε r k B T between 4 and 1/7 (in units of σ). The variation of the temperature allowed covering the wide range of chain states from strongly swollen chains at high temperature to strongly collapsed chains at low temperature. For each temperature, the systems were equilibrated for about 10 6 time steps and then the production runs were performed. Typically, each production run was from 10 6 to 10 7 time steps, depending on the temperature. Throughout the simulation, snapshots of the systems were collected; then using these snapshots, we generated movies illustrating the structural features of the simulated systems.
In the present study we focus on equilibrium properties. They include the mean-square radius of gyration of the chains, the time-averaged potential energies, and the partial pair correlation functions.

Single-chain conformation and coil-to-globule transition
First of all we estimate the reference temperature at which the coil-to-globule transition is observed for the simulated systems. Figure 1 shows the time-averaged squared radius of gyration R 2 g of 128-unit polyions as a function of temperature for z = 1, 2, and 4. An apparent transition temperature, T * , can be associated with the inflection point on the R 2 g (T ) curve. The derivatives ∂ R 2 g (T ) /∂T exhibit pronounced peaks at T * . After a least-squares fitting of the simulation data and subsequent differentiation, we found T * = 1.6, 2.1, and 3.5 at z = 1, 2, and 4, respectively,  for energy parameters ε αβ used in the simulation. We emphasize that the critical temperatures T * estimated here should not necessarily be taken as asymptotic ones (N → ∞); moreover, they should depend upon the concentration of polymer in the solution. Nevertheless, the value of T * should be relevant for the purposes of our study as a reference temperature. At T > T * , the charged copolymer chains are in an expanded coil state and collapse when temperature becomes lower than T * , indicating the formation of dense globules. From the data presented above, we conclude that the transition point shifts towards higher temperature with increasing counterion charge. Also, as seen from figure 1, in the low-temperature region, T < 2.0, the value of R 2 g increases with increasing z while for larger temperatures we find an opposite tendency. This observation will be explained below. We will see that such a behaviour is directly connected with the formation of interchain aggregates of different size, depending on the counterion charge and temperature.

Thermodynamic properties
During the simulations, the average potential energy of the systems was calculated. Figure 2 shows the temperature dependences of the reduced excess energy per particle, U = U/εN (where U is the total potential energy calculated at each time-step and averaged over the full simulation) for different z. For all the systems considered here, the reduced energy U is strongly negative, corresponding to a poor solvent condition, and becomes more negative as the temperature is decreased and the counterion charge increased. Figure 3 demonstrates the temperature dependences of the reduced Coulomb energy U C calculated for different z. As seen, the behaviour of U C is more complex compared to U. For the system containing monovalent counterions, the Coulomb energy U C initially increases with decreasing T , achieves a maximum value at T ≈ 1.5, and then begins to drop, showing a nonmonotonic behaviour in the low-temperature region. What is relevant here is the fact that for all  the temperatures considered in the present study the Coulomb energy remains positive at z = 1. Therefore, the negative values of U seen in figure 2 for this system are due to strong short-range attractive interactions between chain segments. In the presence of di-and tetravalent counterions, we observe a more smooth dependence of U C on T with negative regions at low temperature. In order to explain this behaviour, we discuss the total Coulomb energy as a sum of the following three terms: where U mm is the energy of intra-and interchain monomer-monomer interactions, U cc is the energy of counterion-counterion interactions, and U mc denotes the (negative) energy of polyioncounterion cross interactions. In the high-temperature region, 5.0 < T < 7.0, polyelectrolyte chains have expanded conformations (see figure 1) and counterions are distributed in the system almost uniformly. In this case, the U cc and U mc terms change weakly and the value of U C remains almost constant. As the temperature is decreased, we observe an increase in U C at z = 1 in the temperature region 1.5 < T < 3.0, while in the presence of multivalent counterions the value of U C decreases monotonically. In this temperature region, the behaviour of each individual polyion is mainly governed by the non-ionic short-range attraction between chain beads, resulting in the formation of more compact chain conformations. It is clear that in this case, the distances between charges on each individual chain become shorter, thereby causing an increase in the U mm term. Such a chain contraction results in the growth of the overall electrostatic energy for the z = 1 system and in a more slow decrease in U C in the presence of divalent counterions. For the system containing tetravalent counterions, the negative U mc term dominates. As has been mentioned in section 1, at z = 1, the individual proteinlike globules can form interglobular aggregates in which each globule maintains its morphological integrity and no large-scale aggregation is observed. Due to the fact that the distances between non-compensated charges on single collapsed chains become shorter, this process is accompanied by a jump-like rise in U mm at T ≈ 0.45 thus causing a weak increase in U C . However, such non-monotonic behaviour is not observed when z = 2 or 4. It means that in the presence of multivalent counterions the processes of globule formation and aggregation proceed in a way different from that for the system with monovalent counterions.

Counterion condensation
In order to estimate the critical temperature of counterion condensation, T c , we can use the radius R * , as a measure for the reduced charge of the chain [8,11]. This quantity is defined as the average radius of a sphere around a monomer unit that contains one counterion. The value of R * can be calculated from the equation where ρ c is the number density of counterions and g mc (r) defines the probability of finding a counterion within a spherical shell of volume 4πr 2 dr at a distance r from a given monomer unit, relative to that expected from their random distribution in the system. Figure 4 presents R * as a function of temperature for three values of z. As expected, the value of R * decreases with decreasing temperature and tends to approach unity in the T → 0 limit. The derivative ∂R * /∂T has a rather sharp peak at T c ≈ 0.45 (l B ≈ 2.22) for z = 1, at T c ≈ 1.5 (l B ≈ 0.67) for z = 2 and at T c ≈ 3.0 (l B ≈ 0.33) for z = 4. These temperatures can be considered as the critical temperatures of counterion condensation for the corresponding values of z. From the data shown in figure 1 we can conclude that as the temperature is progressively decreased, the attractive intrachain interaction is so strong that the copolymer chains collapse before counterion condensation starts, i.e., T c < T * for all the z values considered here. Note that for the z = 1 system, the estimated critical temperature T c is considerably lower as compared to the one predicted by the Manning theory [32]. The difference may be due to the fact that in our system counterions mainly condense on the surface of dense globules, which have nearly spherical shape [8]. It is also interesting to note that in the PRISM study of Zherenkova et al [11] the temperature of counterion condensation for charged proteinlike copolymers was found to be higher than for charged random-block copolymers with the same average block length thus indicating the influence of the primary structure on counterion condensation.

Correlation functions
During the simulations, the local structure of the systems was examined via the following paircorrelation functions (PCF): the monomer-monomer PCF, g mm (r); the monomer-counterion PCF, g mc (r); the counterion-counterion PCF, g cc (r); and the polyion-polyion PCF, g pp (r). Note  that g pp (r) describes the spatial distribution of centres-of-mass of polyions. Below, we will focus mainly on the features of g mc (r), g cc (r), and g pp (r). Figures 5 and 6 present the functions g mc (r) and g cc (r) calculated for z = 1, 2, and 4 at a few different temperatures. We note that the first maximum of the monomer-counterion and counterion-counterion PCFs becomes more pronounced with decreasing temperature. The temperature changes observed for g mc (r) show that the local density of counterions around polyelectrolyte chains grows, indicating counterion condensation (see also figure 4). It is seen from figure 5 that in the region of small r (r R 2 g 1/2 ), the first peak of g mc (r) becomes more pronounced as the counterion charge is increased. It means that larger and more compact interchain aggregates are formed in the system containing multivalent counterions. The main maximum of g cc (r) becomes also more pronounced with increasing z and shifts towards larger distances, thereby reflecting the change in the degree of counterion condensation (figure 6).
Using the g pp (r) function, we can estimate how close two charged copolymer chains may come to each other. As an example, figure 7 shows the function g pp (r) for a few different temperatures at z = 1 and 4. Although the averaging was performed over very long runs, there is still substantial statistical noise in the simulation results, especially at short distances and low temperatures. For freely interpenetrating chains, g pp (r) is expected to be constant and equal to unity for all chain-chain distances. However, we see that it exhibits non-random structure even at high temperatures when the chains are in coil-like state. For the high-temperature regime, we observe the development of a correlation hole with g pp (r) = 0 for chain-chain distances r smaller than roughly 10σ. For the z = 1 system, this distance remains practically the same down to T = 0.5 when the polyions are in a dense globular state (see figure 7(a)). This fact indicates that, in the presence of monovalent counterions, charged proteinlike chains having a globular conformation cannot come closer together than this distance. Indeed, in this case the potential of mean force defined as v pp (r) = −k B T ln(g pp (r)) becomes practically infinite. Thus, we can draw the following conclusion: under conditions corresponding to a poor solvent, when the chains are in a globular state, there is a very strong repulsive interaction pushing proteinlike globules away from each other and thus preventing their merging into a single large cluster. The same conclusion has been drawn in [8]. For systems with multivalent counterions, the situation changes dramatically. In the low-temperature region, the correlation hole disappears, and the main maximum of g pp (r) moves to shorter distances with decreasing T . Clearly, these features of the intermolecular pair correlation function are directly connected with a different three-dimensional organization of the collapsed polyions in a poor solvent. Figure 8 presents the position of the main maximum r (1) pp of g pp (r) as a function of temperature for three values of z. To accurately determine the value of r (1) pp , we used the best fit of the g pp (r) functions in the range of their main peaks, all of which were chosen to be Lorentzian in shape. We find that the r (1) pp value shows a sharp decrease in the low-temperature region. It should be noted that the shortest chain-chain distance r (1) pp ≈ 5.4σ observed at z = 1 is larger than the characteristic size of the polyion at the lowest temperature considered here ( R 2 g = 2.31 at T = 0.25) while the analogous shortest distances r (1) pp found for the z = 2 and z = 4 systems are smaller than the corresponding characteristic chain sizes ( R 2 g = 2.54 and 3.24, respectively; see figure 1). Such a behaviour of r (1) pp can be explained as follows. In the presence of monovalent counterions, the distance between the centres-of-mass of the individual globules forming aggregates is larger than the mean chain size due to the fact that in this case, neighbouring globules, in principle, can touch each other or even can stick together, but they cannot interpenetrate [8]. At the same time, in the systems containing multivalent counterions, the chain aggregation proceeds in such a way that several charged copolymers entangle with each other. As a result, the distance between the centres-of-mass of the polyions entering the same multichain aggregate turns out to be smaller than the mean chain size.

Chain aggregation
The aggregation of macromolecules into intermolecular clusters can be described, to first order, using the general framework of the reversible coagulation process [33]. We consider a volume V , which contains n macromolecules. At any time, some fraction of the n chains can be gathered into interchain aggregates of different sizes so that in the system there are n 1 free chains with concentration [A 1 ], n 2 dimers with concentration [A 2 ], and so on. They must satisfy the condition: i n i = n. The aggregation in a dilute solution can be represented by the multipleequilibrium model [34]- [37] described by Association and dissociation is an equilibrium process leading to a size distribution of the intermolecular aggregates. To describe this process, we analyse the average aggregation numbers. The average number m of polyions entering intermolecular aggregates can be estimated from the centre-of-mass pair correlation function, g pp (r). The average number of chains n(r) around a given chain within a distance r is defined as where ρ chain is the number density of the chains, n/V . We assume that two chains belong to an interchain aggregate if their centres-of-mass are located closer to each other than some characteristic 'aggregation distance' r a . In other words, 1 + n(r a ) chains form a connected cluster so that m = 1 + n(r a ). A polyion not belonging to an aggregate is designated as 'unimer'. Certainly, the choice of r a for a multichain system is to some extent arbitrary. For simplicity, we have chosen the parameter r a to be half the characteristic distance between noninteracting particles uniformly distributed at the number density of n/L 3 box . In our case, this gives r a = 0.5L box / 3 √ n = 7.8σ. The average number of copolymer chains per aggregate, m, found in this way is shown in figure 9 as a function of temperature for systems containing counterions of different valency. Obviously, without any attraction between chain segments the time-averaged cluster size is determined by random contacts between the chains and depend primarily on their number density. In this case, one can expect that m ≈ 1 for a dilute solution of non-aggregated chains. The formation of intermolecular aggregates would result in m > 1.
The data presented in figure 9 demonstrate that, regardless of the valency of the counterion, the average aggregate size increases as the temperature is decreased. This is, of course, a quite expected result. One can conclude that for all the systems, the aggregation process becomes well pronounced in the temperature region T T c , where T c is the critical temperature of counterion condensation estimated above. We emphasize that in this case, the chains adopt compact conformations (see figure 1). Although the existence of separated single-chain globules is a prevalent structural motif for the z = 1 system, nevertheless, we observe the occasional formation of intermolecular aggregates: at T = 0.25, the globules ('unimers') can stick together, forming pairs and (very seldom) trimers and larger aggregates. In the presence of di-and tetravalent ions, the average aggregate size is considerably larger. Nevertheless, for all the systems studied here the interchain aggregates have a well-defined finite size even at the lowest temperature; that is, no large-scale aggregation is observed.

Internal aggregate structure
As has been discussed elsewhere [8], in the case of monomolecular proteinlike globules formed at low temperature in the presence of monovalent cations, there is a quite well-defined inner (central) region of hydrophobic chain segments surrounded by an outer layer filled with charged segments and cations. Thus, we do not observe a homogeneous structure of the globules with chain groups homogeneously distributed inside the globule. One may say that an intramolecular microphase separation of H and P chain groups takes place inside a proteinlike globule.
To clarify the internal equilibrium structure of the interchain aggregates formed in the presence of multivalent counterions, we calculated the density profiles for H and P polymer segments and for counterions. As an example, we show in figure 10 the normalized radial distribution functions W(r) calculated for z = 4 and T = 1, where the centre of mass of the aggregate is taken as the origin. As can be seen, there is no sharp interface between the aggregate core and counterions. Moreover, we observe penetration of the counterions into the aggregate core, although they show somewhat inhomogeneous distribution, tending to localize in outer and inner aggregate regions. Partly, this is due to the size and/or shape fluctuations of the aggregate as a whole. It should be stressed that even at low temperatures, there is no exact neutralization of negative and positive charges in the volume occupied by multichain aggregates. A certain fraction of counterions always remains in solution. Therefore, there is redistribution of counterions between the interior of the aggregates and the solution outside. Also, we note that charged chain segments tend to be on the surface of the aggregates, thereby screening their hydrophobic interior from the solvent. Figure 11 illustrates a typical instantaneous picture (snapshot) of an intermolecular aggregate formed at T = 0.25 and z = 4. The aggregate is built up from 7 polyions, which are entangled together, and 110 tetravalent counterions so that the net charge of the aggregate is Q = −8. For visual clarity, one of the chains is depicted by thick rods. It is seen that tetravalent counterions penetrate into the interior of the nearly spherical aggregate and also condense on its surface. The aggregate surface is covered with polar chain sections. Some fraction of counterions is floating in the immediate vicinity of the polyions. All these features are also seen clearly in the computer-generated movie. Computer-generated movie illustrating the structural features of the aggregate and the distribution of counterions. The movie is generated from subsequent snapshot pictures of a MD simulation. We used the PovRay software (http://www.povray.com) for rendering frames and TMPGEnc software (http://www.tmpgenc.net) to convert the rendered frames into a movie in mpeg1 format. To play the movie, a DivX codec (http://www.divx.com) may have to be installed.

Discussion
We begin with a short discussion of the mechanism responsible for the formation of finite-size globular aggregates and their stabilization under poor solvent conditions. In our model, there are two driving forces causing the formation of compact conformations of charged copolymer chains: the short-range (van der Waals type) attraction and the counterionmediated attraction between charged chain segments. In recent years, the origin of counterionmediated attraction has been extensively discussed in the literature (see, e.g. [38]- [44]). About 10 years ago, Ray and Manning [45] suggested that two like-charged objects (e.g., rigid rods) can share condensed counterions in such a way to allow them to form what is analogous to a chemical bond. This 'bridging type' model has been extensively used to explain the counterioninduced precipitation of polyelectrolytes [12,46,47] and other related phenomena. When the temperature is reduced, we observe for our model the counterion condensation (figure 4) that makes the counterion-mediated attraction stronger. As we discussed earlier, this effect is more pronounced for counterions with higher charge.
The stabilization of the finite-size aggregates (both single-chain and multichain) is due to the interplay between the attractive forces and the long-range Coulomb interactions together with the contribution from the translational entropy of mobile counterions. The point is that only a fraction of the charges on the polyions are neutralized by condensed counterions because the temperature is not zero, so that each copolymer chain still carries a net negative charge. Therefore, when the formation of the aggregates occurs, their net charge (i.e., the sum of the net charges of the individual chains in the aggregate) turns out to be non-zero for all aggregate sizes. The electrostatic repulsion due to the net-charge, whose range is set by the screening length under given conditions, is the first factor responsible for the stabilization of finite size aggregates. Obviously, the repulsion is longer in range than the attraction. Also, our simulations show that for the HP copolymers having proteinlike primary structure, the significant fraction of charged chain segments is concentrated on the aggregate surface (see figures 5, 10 and 11). Such a morphology assists the stabilization. On the other hand, it seems evident that this will not be the case for polyions with a purely random primary structure or for regular copolymers with an alternating distribution of charged and neutral groups along the chain. In the counterion condensation regime, the remaining (non-condensed) counterions are redistributed between the outer solution and the volume occupied by the polyions included into aggregates. The inhomogeneous distribution of counterions results in translational entropy losses. This is the second ingredient of the stabilization mechanism. For small clusters, the long-range Coulomb repulsion is more important, while for large interchain aggregates the stabilization should originate mainly from the counterion entropy decrease. Indeed, as the aggregates grow and occupy a larger fraction of the available volume in the system, a larger fraction of the free counterions will be found in the aggregates, thereby leading to an increase in the corresponding free energy connected with the translational entropy losses. Also, it should be noted that the counterion mobility inside an aggregate is reduced when the chains contract significantly and the aggregates become denser. All these arguments explain the existence of the finite-size aggregates built up from charged copolymers in a poor solvent.
The stabilization mechanism described above is practically the same as that suggested in [21]- [24]. In particular, the theory developed in this work predicts that in the dilute solution of associating polyelectrolytes, finite-size clusters having some optimum dimension should exist. Moreover, this stabilization mechanism is universal and should be valid for all attracting polyelectrolytes independent of the nature of the attraction between chain segments and the internal structure of the aggregates. To illustrate this assertion, the authors of [22] present the following qualitative arguments. When the polyions aggregate into a multichain cluster, the intermolecular potential energy of attraction per chain is a monotonically decreasing function of the aggregation number, m. It is clear that this energy should decrease up to some finite negative constant value for the macroscopic cluster in the m → ∞ limit. On the other hand, the sum of the electrostatic repulsive energy and the counterion entropy decrease (per chain) is a monotonically increasing function of m. As a result, there should be a minimum of the total free energy at some finite value of m that corresponds to the formation of optimum finite-size aggregates for the given conditions. An analogous stabilization mechanism was also predicted by Borue and Erukhimovich [48] and then studied in detail by Joanny and Leibler [49].
We can add a few words about the behaviour of the chain size observed in figure 1 for the low-temperature region. Let us compare the conformation of a polyion entering a single-chain and multichain aggregate. In the single-chain aggregate (i.e., in a dense polymer globule), every strongly collapsed polyion has a size of the order of N 1/3 . In the large multichain aggregate, however, each flexible-chain macromolecule is entangled with other chains, which form a surrounding essentially similar to that characteristic of a polymer melt. In the m → ∞ limit, each flexible chain should approximately behave as a Gaussian chain with the size of the order of N 1/2 [50], i.e., its size increases. That is why the chains included into larger aggregates formed in the presence of multivalent ions have larger sizes than those observed for small single-chain aggregates (see figure 1). The chain expansion is also clearly seen in figure 11.

Conclusion
We have presented the results of molecular dynamics simulations of dilute salt-free polyelectrolyte solutions containing charged proteinlike copolymers and multivalent counterions. In this study, the main attention has been paid to the influence of differently charged counterions on the aggregation of the copolymers in a poor solvent. It was found that the aggregation occurs at higher temperature as the counterion charge is increased. In the presence of multivalent counterions, the coil-to-globule transition and counterion condensation also shift towards higher temperatures. Although condensed counterions of higher valency induce stronger interchain attraction and weaker repulsion, the aggregation of proteinlike copolymers is not accompanied by large-scale aggregation, that is, no macroscopic phase separation is observed. Instead, we have observed the formation of thermodynamically stable finite-size aggregates built up from several entangled copolymer chains. This stabilization mechanism can be explained by the interplay between the short-range attractive forces and the long-range Coulomb interactions together with the contribution from translational entropy of mobile counterions. ∞ x e −t 2 dt, and α is an inverse length (the so-called Ewald free parameter). The real space and Fourier space forces on the charged particle i, obtained by differentiating the potential energies defined by (A.2) and (A.3), are given by √ π e −α 2 |r ij +mL box | 2 + erfc(α|r ij + mL box |) |r ij + mL box | r ij + mL box |r ij + mL box | 2 , (A.5) The self-energy term (A.4) does not depend on r and k and hence does not produce forces on particles.
To take into account long-range electrostatic interactions between charges, we used the highly optimized particle mesh Ewald (PME) method [52,53] and the fast Fourier transform algorithm. The truncation distance for electrostatic interactions in equation (A.2) was set to R C = 10σ, and long-range corrections to the potential energy were made. The Ewald convergence parameter was α = 0.3766σ −1 . The order of B-splines was n B = 5 for all systems. The number of grid points in each direction was 72. These parameters ensured consistency between the Coulomb energy obtained through the pair potential and the Coulomb energy calculated through the virial.