A Excluded volume contributions

Self-assembly of short DNA duplexes: from a coarse-grained model to experiments through a theoretical link


Short blunt-ended DNA duplexes comprising 6 to 20 base pairs self-assemble into polydisperse semi-flexible chains due to hydrophobic stacking interactions between terminal base pairs. Above a critical concentration, which depends on temperature and duplex length, such chains order into liquid crystal phases. Here, we investigate the self-assembly of such double-helical duplexes with a combined numerical and theoretical approach. We simulate the bulk system employing the coarse-grained DNA model recently proposed by Ouldridge et al. [ J. Chem. Phys. 134, 08501 (2011) ]. Then we evaluate the input quantities for the theoretical framework directly from the DNA model. The resulting parameter-free theoretical predictions provide an accurate description of the simulation results in the isotropic phase. In addition, the theoretical isotropic-nematic phase boundaries are in line with experimental findings, providing a route to estimate the stacking free energy.

I Introduction

Self-assembly is the spontaneous formation through free energy minimization of reversible aggregates of basic building blocks. The size of the aggregating units, e.g. simple molecules, macromolecules or colloidal particles, can vary from a few angströms to microns, thus making self-assembly ubiquitous in nature and of interest in several fields, including material science, soft matter and biophysics (1); (2); (3); (4); (5). Through self-assembly it is possible to design new materials whose physical properties are controlled by tuning the interactions of the individual building blocks (6); (7); (8); (9).

A relevant self-assembly process is the formation of filamentous aggregates (i.e. linear chains) induced by the anisotropy of attractive interactions. Examples are provided by micellar systems (10); (11); (12), formation of fibers and fibrils (13); (14); (15); (16), solutions of long duplex B-form DNA composed of to base pairs  (17); (18); (19); (20), filamentous viruses (21); (22); (23); (24); (25), chromonic liquid crystals (26) as well as inorganic nanoparticles (27).

If linear aggregates possess sufficient rigidity, the system may exhibit liquid crystal (LC) phases (e.g. nematic or columnar) above a critical concentration. In the present study we focus on the self-assembly of short (i.e. 6 to 20 base pairs) DNA duplexes (DNADs) (28); (29); (30) in which coaxial stacking interactions between the blunt ends of the DNADs favor their aggregation into weakly bonded chains. Such a reversible physical polymerization is enough to promote the mutual alignment of these chains and the formation of macroscopically orientationally ordered nematic LC phases. At present, stacking is understood in terms of hydrophobic forces acting between the flat hydrocarbon surfaces provided by the paired nucleobases at the duplex terminals, although the debate on the physical origin of such interaction is still active and lively (31); (32). In this respect, the self-assembly of DNA duplexes provides a suitable way to access and quantify hydrophobic coaxial stacking interactions.

In order to extract quantitative informations from DNA-DNA coaxial stacking experiments, reliable computational models and theoretical frameworks are needed. Recent theoretical approaches have focused on the isotropic-nematic (I-N) transition in self-assembling systems  (33); (34), building on previous work on rigid and semi-flexible polymers (35); (36); (37); (38); (39); (40); (41); (42); (43). In a recent publication (44) we investigated the reversible physical polymerization and collective ordering of DNA duplexes by modeling them as super-quadrics with quasi-cylindrical shape (45) with two reactive sites (46); (47) on their bases. Then we presented a theoretical framework, built on Wertheim (48); (49); (50) and Onsager (51) theories, which is able to properly account for the association process.

Here, we employ this theoretical framework to study the physical properties of a realistic coarse-grained model of DNA recently proposed by Ouldridge et al. (52), where nucleotides are modeled as rigid bodies interacting with site-site potentials. Following Ref. (44), we compute the inputs required by the theory, i.e. the stacking free energy and the DNAD excluded volume, for the Ouldridge et al. model (52). Subsequently we predict the polymerization extent in the isotropic phase as well as the isotropic-nematic phase boundaries.

To validate the theoretical predictions, we perform large-scale molecular dynamics (MD) simulations in the NVT ensemble of a bulk system comprising nucleotides, a study made possible by the computational power of modern Graphical Processing Units (GPUs). The parameter-free theoretical predictions provide an accurate description of the simulation results in the isotropic phase, supporting the theoretical approach and its application in the comparison with experimental results.

The article is organized as follows. Section II provides details of the coarse-grained model of DNADs, of the MD computer simulations and presents a summary of the theory. Section IV describes the protocols implemented to evaluate the input parameters required by the theory via MC integrations for two DNADs. We also discuss some geometrical properties of the bonded dimer configurations. We then compare the theoretical predictions with simulation and experimental results. Finally, in Section V we discuss estimates for the stacking free energy and present our conclusions.

Ii Methods

ii.1 Model

Figure 1: (a) Schematic representation of the coarse graining of the model for a single nucleotide. (b) Model interaction sites. For the sake of clarity, stacking and hydrogen-bonding sites are highlighted on one nucleotide and the base repulsion site on the other. For visualization reasons, in the following strands will be shown as ribbons and bases as extended plates as depicted in (c). (d) A base-pairs DNAD in the minimum energy configuration. (e) An equilibrium configuration of the same object at . The nucleotides at the bottom are not bonded, the so-called fraying effect. (f) A chain of length extracted from a simulation at .

We implement a coarse-grained model for DNA recently developed by Ouldridge et al. (53); (52). In such model, designed via a top-down approach, each nucleotide is described as a rigid body (see Figure 1). The interaction forms and parameters are chosen so as to reproduce structural and thermodynamic properties of both single- (ssDNA) and double- (dsDNA) stranded molecules of DNA in B-form. All interactions between nucleotides are pairwise and, in the last version of the model (52), continuous and differentiable. Such feature makes the model convenient foe MD simulations.

The interactions between nucleotides account for excluded volume, backbone connectivity, Watson-Crick hydrogen bonding, stacking, cross-stacking and coaxial-stacking. The interaction parameters have been adjusted in order to be consistent with experimental data (52); (54); (55). In particular, the stacking interaction strength and stiffness have been chosen to be consistent with the experimental results reported for 14-base oligomers by Holbrook et al. (55). Hydrogen-bonding and cross-stacking potentials were adjusted to give duplex and hairpin formation thermodynamics consistent with the SantaLucia parameterization of the nearest-neighbor model (54). Interaction stiffnesses were also further adjusted in order to provide correct mechanical properties of the model, such as the persistence length. The model does not have any sequence dependence apart from the Watson-Crick pairing, meaning that the strength of the interactions acting between nucleotides is to be considered as an average value. In addition, the model assumes conditions of high salt molarity (). In this model, the coaxial-stacking interaction acts between any two non-bonded nucleotides and is responsible for the duplex-duplex bonding. It has been parametrized (56) to reproduce experimental data which quantify the stacking interactions by observing the difference in the relative mobility of a double strand where one of the two strands has been nicked with respect to intact DNA (57); (58) and by analyzing the melting temperatures of short duplexes adjacent to hairpins (59).

To cope with the complexity of the model and the large number of nucleotides involved in bulk simulations, we employ a modified version of the CPU-GPU code used in a previous work (60), and extend it to support the force-fields (61). Harvesting the power of modern Graphical Processing Units (GPUs) results in a -fold speed-up. The CPU version of the code, as well as the Python library written to simplify generation of initial configurations and post-processing analysis, will soon be released as free software (62).

ii.2 Bulk simulations

To compare numerical results with theory, we perform Brownian dynamics simulations of dsDNA molecules made up by nucleotides each, i.e. cylinder-like objects with an aspect ratio of (see Figure 1 (d)). The integration time step has been chosen to be in simulation units which corresponds, if rescaled with the units of length, mass and energy used in the model, to approximately seconds.

We study systems at three different temperatures, namely , and , and for different concentrations, ranging from to . The state point, despite being far from the experimentally accessed , is here investigated to test the theory in a region of the phase diagram where the degree of association is significant. To quantify the aggregation process we define two DNADs as bonded if their pair interaction energy is negative. Depending on temperature and concentration, we use MD steps for equilibration and MD steps for data generation on NVIDIA Tesla C2050 GPUs, equivalent to .

Figure 2: Snapshots taken from simulations at . At low concentrations (, top) chain formation is negligible and the average chain length is approximately . As the concentration is increased (, bottom), DNADs start to self-assemble into chains and the average chain length increases.

Iii Theory

We build on the theoretical framework previously developed to account for the linear aggregation and collective ordering of quasi-cylindrical particles (44). Here, we provide a discussion of how such a theory can be used to describe the reversible chaining and ordering of oligomeric DNADs at the level of detail adopted by the present model. According to Ref. (44), the free energy of a system of equilibrium polymers can be written as


where is the volume of the system, ( is the number density of monomers) is the packing fraction, is the number density of chains of length , normalized such that , is the volume of a monomer, , as discussed in the subsection IV.2, is a parameter which depends on the free energy associated to a single bond and is the excluded volume of two chains of length and . is the Parsons-Lee factor (63)


and (43) accounts for the orientational entropy that a chain of length loses in the nematic phase (including possible contribution due to its flexibility). The explicit form for can be found in Ref.  (44).

The free energy functional (Eq. 1) explicitly accounts for the polydispersity inherent in the equilibrium polymerization using a discrete chain length distribution and for the entropic and energetic contributions of each single bond through the parameter .

Isotropic phase

In the isotropic phase and the excluded volume can be written as follows (see Appendix A):


where the parameters and can be estimated via MC integrals of a system composed by only two monomers (see Appendix A) and is the aspect ratio of the monomers. We assume that the chain length distribution is exponential (44) with an average chain length




With this choice for the free energy in Eq. (1) becomes:


Minimization of the free energy in Eq. (6) with respect to provides the following expression for the average chain length :


Nematic phase

In the nematic phase the monomer orientational distribution function depends explicitly on the angle between the particle and the nematic axis, i.e. the direction of average orientation of the DNAD, since the system is supposed to have azimuthal symmetry around such axis. We assume the form proposed by Onsager (51), i.e.:


where controls the width of the angular distribution. Its equilibrium value is obtained by minimizing the free energy with respect to . As discussed in Appendix A, we assume the following form for the excluded volume in the nematic phase:


where the term is the end-midsection contribution to the excluded volume of two hard cylinders (see Appendix B) and


In Eq. (10), is the diameter of the monomer and with are three parameters that we chose in order to reproduce the excluded volume calculated from MC calculations as discussed in Appendix A.

Inserting Eqs. (9) and (4) into Eq. (1) and assuming once more an exponential distribution for one obtains after some algebra:


where . The explicit calculation of the parameters and is explained in Appendices  A and B.

Assuming that the orientational entropy can be approximated with the expression valid for long chains (43), minimization with respect to results in


while using the approximated expression for short chains (43), one obtains


The equilibrium value of is thus determined by further minimizing the nematic free energy in Eq. (11), which has become only a function of . The parameter is related to the degree of orientational ordering in the nematic phase as expressed by the nematic order parameter as follows:


Further refinements of the theory could be obtained by including a more accurate description of the orientational distribution in the proximity of the I-N phase transition, along the lines of Eqs. (40)-(42) of Ref. (44). For the sake of simplicity we have just presented the basic theoretical treatment. However, in the theoretical calculations in Sec. IV we will make use of the refined and more accurate free energy proposed in Ref. (44).

Phase Coexistence

The phase boundaries, at which the aggregates of DNAD are sufficiently long to induce macroscopic orientational ordering, are characterized by coexisting isotropic and nematic phases in which the volume fraction of DNADs are, respectively, and . The number densities and can be calculated by minimizing Eq. (6) with respect to and by minimizing Eq. (11) with respect to and . In addition, the two phases must be at equal pressure, i.e. , and chemical potential, i.e. . These conditions yield the following set of equations:


Iv Results and Discussion

iv.1 Properties of the model

To characterize structural and geometrical properties of the simulated DNADs monomers and aggregates we analyze conformations of duplexes extracted from large-scale GPU simulations (see Figure 2 for some snapshots).

In the following, the volume occupied by a single DNAD of length and double helix diameter () will be calculated as the volume of a cylinder with the same length and diameter, i.e. . When comparing numerical and experimental results with theoretical predictions we use the number of nucleotides in place of () and the concentration instead of the packing fraction , which can be related to the former via:


where is the average mass of a nucleotide. Hence, in the following and will be used in place of and .

First we calculate the dimensions (height and width ) of the DNADs for different and . We observe no concentration dependence on both quantities, while the variation in is negligible (of the order of between DNADs of samples at and ). The effect of this small change does not affect substantially the value of the aspect ratio, which we consider constant () throughout this work.

The geometrical properties of end-to-end bonded duplexes are not well-known since there are no experimental ways to probe such structures. In a very recent work, the interaction between duplex terminal base–pairs has been analyzed by means of large-scale full-atom simulations by Maffeo et al. (64). They found that blunt-ended duplexes (i.e. duplexes without dangling ends) have preferential binding conformations with different values of the azimuthal angle , defined as the angle between the projections onto the plane orthogonal to the axis of the double helix of the vectors connecting the O5’ and O3’ terminal base pairs. They report two preferential values for , namely and .

Figure 3: Probability distributions for (a) the azimuthal angle and (b) the end-to-end distance .

In the present model the continuity of the helix under end-to-end interactions is intrinsic in the model and the azimuthal angle probability distribution is peaked around a single value (see Figure 3(a)). This is very close to the theoretical value given by the pitch of the B-DNA double helix. The qualitative difference between the conformations of bonded DNADs found in this work and in Ref.  (64) should be addressed in future studies describing the coaxial end-to-end interaction in a more proper way.

In addition, we calculate the average distance between the centres of masses of the terminal base pairs. Figure 3(b) shows , the probability distribution of . is peaked at , whereas Maffeo et al. (64) found an average distance of . This difference can be understood in terms of the effect of the salt concentration which, being five times higher than the one used in Ref. (64), increases the electrostatic screening, thus effectively lowering the repulsion between DNA strands.

The effect of the temperature is small, as lowering leads only to more peaked distributions for both and (and a very small shift towards smaller angles for ) but does not change the overall behavior.

iv.2 Stacking free energy and excluded volume

In this section we discuss the procedure employed to evaluate the input quantities required by the theory, namely and . To this aim we perform a Monte Carlo integration over the degrees of freedom of two duplexes. is defined as (44)


where (65)


Here is the vector joining the center of masses of particles and , is the orientation of particle and represents an average taken over all the possible orientations. is the bonding volume, defined here as the set of points where the interaction energy between duplex and duplex is less than . To numerically evaluate we perform a MC integration using the following scheme:

  1. Produce an ensemble of 500 equilibrium configurations of a single duplex at temperature .

  2. Set the counter and the energy factor .

  3. Choose randomly two configurations and from the generated ensemble.

  4. Insert a randomly oriented duplex in a random position in a cubic box of volume . Insert a second duplex in a random position and with a random orientation. Compute the interaction energy between the two duplexes and and, if , update the energy factor, . Increment .

  5. Repeat from step 3, until converges within a few per cent precision.

The employed procedure to compute is fairly similar except that it is performed for duplexes with a various number of bases (i.e. with different ) and the quantity counts how many trials originate a pair configuration with (i.e. in step , ). In the nematic case, the orientations of the duplexes are extracted randomly from the Onsager distribution given by Eq. 8. With such procedure,


We calculate for values of , ranging from to (see Appendix B). Since the and dependences of Eqs. (3) and (9) are the same and the dependence of the numerically calculated on the shape of DNADs is negligible, the evaluation of the excluded volume as a function of provides the same information as the evaluation of as a function of .

Fig. 4 shows for all investigated in a vs plot. A linear dependence properly describes the data at the three . An alternative way to evaluate is provided by the limit of Eq. (7). Indeed in the low density limit and are related via the following relation:


Therefore it is also possible to estimate by extrapolating the low density data for at , and . The results, also shown in Fig. 4, are in line with the ones obtained through MC calculations. The Arrhenius behavior of suggests that bonding entropy and stacking energy are in first approximation independent. The coaxial stacking free energy is related to as follows


Substituting the fit expression provided in Fig. 4 for results in a stacking free energy at a standard concentration of DNADs and comprising a bonding entropy of and a bonding energy of .

Figure 4: calculated with the procedures described in Sec. IV.2 for all investigated .

iv.3 Isotropic phase: comparing simulation results with theoretical predictions

Fig. 5 shows the concentration dependence of calculated from the MD simulation of the system. The average chain length increases progressively on increasing . The figure also shows the theoretical predictions calculated by minimizing the isotropic free energy in Eq. (6) with respect to using the previously discussed estimates for and . The theoretical results properly describe the MD simulation data up to concentrations around , which corresponds to a volume fraction . In Ref. (44) similar observations have been made and the discrepancy at moderate and high has been attributed to the inaccuracy of the Parsons decoupling approximation. The values calculated using the excluded volume of two hard cylinders (HC) are also reported, to quantify the relevance of the actual shape of the DNA duplex. Indeed the HC predictions appreciably deviate from numerical data beyond mg/ml.

Figure 5: Average chain length in the isotropic phase at low concentration. Symbols are numerical results and dashed lines are theoretical predictions. Dotted lines are theoretical predictions using the excluded volume of HCs (see Appendix B).

iv.4 Phase Coexistence: Theoretical predictions

A numerical evaluation of the phase coexistence between the isotropic and the nematic phases for the coarse-grained model adopted in this study is still impossible to obtain given the current computational power. We thus limit ourselves to the evaluation of the I-N phase coexistence via the theoretical approach discussed in Sec. III. Fig. 6 shows the theoretical phase diagram in the - plane for and . As expected, both and decrease on increasing , since the increase of the number of basis result in a larger aspect ratio. On decreasing , theory predicts a decrease of and a similar decrease of , resulting in an overall shift of the I-N coexistence region toward lower values. This trend is related to the increase of the average chain length with increasing (see Fig. 4).

Fig. 6 also shows the phase boundaries calculated using the excluded volume of two hard cylinders. Assimilating DNADs to hard cylinders results in a widening of the isotropic-nematic coexistence region.

Figure 6: I-N phase diagram in the vs plane for (top) and (bottom). Dotted lines are theoretical phase boundaries calculated using the excluded volume of HCs (see Appendix B).

iv.5 Comparison between theory and experiments

The theoretical predictions concerning the isotropic-nematic coexisting concentrations can be compared to the experimental results reported in Refs. (28) and (30) for blunt-ended DNADs.

Figure 7 compares the experimentally determined nematic concentrations at coexistence with the values calculated from the present model for . Despite all the simplifying assumptions and despite the experimental uncertainty, the results provide a reasonable description of the dependence of . The experimental data refer to different base sequences and different salt concentrations. According to the authors is affected by an error of about . In particular for the case the critical concentrations for distinct sequences shows that blunt-end duplexes of equal length but different sequences can display significantly different transition concentrations. Hence, for each duplex length, we consider the lowest transition concentration among the ones experimentally determined, since this corresponds to the sequence closest to the symmetric monomer in the model. Indeed the dependence of on the DNADs sequence is expected to be larger for the shortest sequences, i.e. , for which DNAD bending could be significant (30). Unfortunately, quantitative experimental data on this bending effect are still lacking. In general it is possible that for (for which a large number of sequences have been studied, see Fig. 7), would be corrected to lower values if a larger number of sequences were explored. For more details on this phenomenon, we refer the reader to the discussions in Refs. (30); (66); (44).

The overestimation of the phase boundaries for with respect to experimental results suggests that the DNA model of Ouldridge et al. (52) overestimates the coaxial stacking free energy. Such discrepancy can perhaps be attributed to the restricted number of microstates allowing for bonding states in the DNA model (52); (56), as discussed in Sec. IV. Indeed, allowing DNADs to form end-to-end bonds with more than one preferred azimuthal angle would increase the entropy of bonding, thus effectively lowering . Allowing for both left- and right-handed binding conformations, a possibility supported by the results of Maffeo et al. (64), would double in Eq. (21) and hence add an entropic contribution equal to to , which would result in a value for . Fig. 7 also shows the theoretical prediction for such upgraded value.

In Fig. 7 theoretical calculations of the I-N transition lines are shown for and at as the upper and lower boundaries of the grey band respectively. To calculate these critical lines we retain the excluded volume calculated in the subsection IV.2 and, given the value of , we evaluate according to Eqs. (17) and (21) for and corresponding to the standard concentration.

The selected points with fall within the grey band shown in Fig. 7 enabling us to provide an indirect estimate of between and . For the points with , where duplex bending might play a role, it would be valuable to have more experimental points corresponding to more straight sequences in order to validate the theoretical predictions.

It is worth observing that for all DNAD lengths the electrostatics interactions are properly screened. For a concentration of NaCl has been added to the solution resulting in a Debye screening length . For all other lengths (i.e. ) we note that at the lowest DNA concentration of corresponding to . Therefore the experimental is always smaller than the excluded volume diameter for the backbone-backbone interaction of our coarse-grained model (52) (), thus enabling us to neglect electrostatic interactions.

On the other hand, if electrostatics interactions are not properly screened the effective aspect ratio for such DNAD sequences would be smaller than the ones used in our theoretical treatment and this would result in a underestimate of . To account for this behavior one should at least have a reasonable estimate of the effective size of DNADs when electrostatics interactions are not fully screened. Moreover, the role of electrostatics interactions can be subtle and not completely accounted for by simply introducing an effective size of DNADs. A possible route to include electrostatics in our treatment can be found in Ref. (20) and it will be addressed in future studies.

Figure 7: Critical nematic concentrations as a function of the number of base pairs per duplex for the present model, calculated theoretically at using the computed stacking free energy (short dashed lines), (long dashed lines), and for experiments (28) (circles and squares). Squares are for different sequences at the same . The grey band has been built considering for an upper bound of and a lower bound of .

iv.6 Comparison with Onsager Theory

The experimental average aggregation numbers are estimated in Refs. (13); (28) by mapping the self-assembled system onto an “equivalent” mono-disperse system of hard rods with an aspect ratio equal to . In Ref. (44) it has been shown that the theoretically estimated isotropic-nematic coexistence lines for the case of polymerizing superquadric particles in the plane, parametrized by the stacking energy, are significantly different from the corresponding Onsager original predictions (as re-evaluated in Ref. (35)). In light of the relevance for interpreting the experimental data, we show in Figure 8 the same curves for the DNA model investigated here. In this model, a clear re-entrant behavior of the transition lines in the plane is observed. The re-entrant behavior occurs for values of the stacking free energy accessed at temperatures between and . We believe that the re-entrancy of the transition lines in the plane is a peculiar mark of the system polydispersity resulting from the reversible self-assembling into chains, and as such it should be also observable experimentally.

Figure 8: Isotropic-nematic coexistence lines in the average aspect ratio and concentration plane for two values of , namely (top) and (bottom). Solid lines indicate theoretical predictions, dashed lines indicate the Onsager original predictions, as re-evaluated in Ref.  (35) for and . Symbols along the isotropic and nematic phase boundaries at coexistence are joined by dotted lines, to indicate the change in concentration and average chain length at the transition.

V Conclusions

In this article, we have provided the first study of a bulk solution of blunt ended DNA duplexes undergoing reversible self-assembly into chains, promoted by stacking interactions. The simulation study, carried out at different concentrations and temperatures, provides a clear characterization of the and dependence of the average polymerization length and an indirect estimate of the stacking free energy. We have provided a theoretical description of the self-assembly process based on a theoretical framework recently developed in Ref. (44). The inputs required by the theory (the DNAD excluded volumes and the stacking free energy) have been numerically calculated for the present DNA model, allowing a parameter free comparison between the molecular dynamic results and the theoretical predictions. Such comparison has been limited to the isotropic phase, due to the difficulties to simulate the dense nematic phase under equilibrium conditions. The description of the isotropic phase is satisfactory: quantitative agreement between theory and simulations is achieved for concentrations up to mg/ml. The stacking free energy value that properly accounts for the polymerization process observed in the molecular dynamics simulations is at a standard concentration 1 M of DNADs and K comprising a bonding entropy of and a bonding energy of .

Theoretical predictions for the I-N transition have been compared with experimental results for several DNA lengths, ranging from 8 to 20 bases. For the model predicts values for which are higher than experimental ones. This suggests that the DNA model employed overestimates . In view of the recent results of Maffeo et al. (64), we speculate that the bonding entropy is underestimated, in agreement with the observation that the probability distribution of the azimuthal angle between two bonded DNADs, which is designed to be single-peaked, is too restraining. In this respect, the present study call for an improvement of the coarse-grained potential (52) in regard to the coaxial stacking interaction.

The value of can also be used as a fitting parameter in the theory for matching with the experimental results, retaining the excluded volume estimates calculated for the coarse-grained DNA model. Such procedure shows that values of the stacking free energy between and are compatible with the experimental location of the I-N transition line. In the work of Maffeo et al., the authors report a quite smaller value of , namely , a value which was confirmed by the same authors by performing an investigation of the aggregation kinetic in a very lengthy all-atom simulation of DNAD with . If such value is selected as input in our theoretical approach (maintaining the same excluded volume term), then one finds mg/ml, a value significantly smaller than the experimental result ( mg/ml). This casts some doubts on the effectiveness of the employed all-atom force-field to properly model coaxial stacking.

Finally, our work draws attention to the errors affecting the estimate of the average chain length via a straightforward comparison of the nematic coexisting concentrations with analytic predictions based on the original Onsager theory for mono-disperse thin rods (28); (13). We have found that such approximation significantly underestimates at the I-N transition concentration . In addition, the theoretical approach predicts a re-entrant behavior of the transition lines in the - plane, a distinct feature of the polydisperse nature of the equilibrium chains.


We thank Thomas Ouldridge, Flavio Romano and Teun Vissers for fruitful discussions. We acknowledge support from ERC-226207-PATCHYCOLLOIDS and ITN-234810-COMPLOIDS as well as from NVIDIA.

Appendix A Excluded volume contributions

Here we further discuss the calculation of the excluded volume term for the present model. Following Ref. (44), the excluded volume is assumed to be the following second order polynomial in and :


where the functions , , describe the angular dependence of the excluded volume. The orientational probability is normalized such that


The three contributions to the excluded volume in Eq. (22) come from end-end, end-midsection and midsection-midsection steric interactions (44) between two chains.

In the isotropic phase the orientational distribution does not have any angular dependence, i.e. , and Eq. (22) reduces to the form


The parameters , and appearing in Eq. (24) can be calculated via MC integration procedures as discussed in the subsection IV.2 and in Ref. (44). We expect that these parameters do not depend on because each DNADs comprises stacked base pairs which are all identical with respect to excluded volume interactions (i.e. they all have the same shape). In particular, the calculated excluded volume of two DNADs is reported in Fig. 9 for different aspect ratios, together with the resulting values for the above parameters.

Figure 9: Excluded volume in the isotropic phase together with analytic approximations. From the linear fit one has and , while we assume

Using the Onsager angular distribution in Eq. (24), the excluded volume in the nematic phase depends also on the parameter , i.e. the general form in Eq. (22) reduces to


Assuming that , and is given by Eq. (10), the three parameters with have to be estimated. For and several values of ( in steps of ) and we calculated numerically the nematic excluded volume for two DNADs. The results are shown in Fig. 10, where we plot vs for various . The dashed lines shown in Fig. 10 are obtained through a two-dimensional fit to numerical data for using Eq. (9) as fitting function.

Appendix B Excluded volume of hard cylinders

Figure 10: Excluded volume as a function of aspect ratio in the nematic phase together with analytic approximations for several . The dashed lines are obtained plotting the function reported in Eq. (9) and setting , and .

For two rigid chains of length and which are composed of hard cylinders (HCs) of diameter and length , can be described by


where , and are the orientations of two HCs and is the complete elliptical integral


The integrals in Eq. (26) can be calculated exactly in the isotropic phase, while in the nematic phase the calculation can be done analytically only for suitable choices of the angular distribution . Here we assume that the angular distribution is given by the Onsager function in Eq. (8).

Using the Onsager orientational function the following approximate expressions for the coefficients , and can be derived (43)




We evaluate numerically the excluded volume in Eq. (26) for many values of and, building on the expressions in Eqs. (28), we perform a fit to this data using the following functions:


The coefficient values resulting from the fitting procedure are , , , , , and .


  1. To whom correspondence should be addressed. Tel: +390649913524; Fax: +39064463158; Email: cristiano.demichele@roma1.infn.it


  1. I. Hamley, Introduction to Soft Matter, Wiley & Sons, 2007.
  2. S. C. Glotzer, Science, 306, 419–420 (2004). Some assembly required.
  3. G. M. Whitesides and M. Boncheva, Proc. Natl. Acad. Sci. USA, 99(8), 4769–4774 (2002). Beyond molecules: Self-assembly of mesoscopic and macroscopic components.
  4. Y.-S. Cho, G.-R.Yi, J.-M. Lim, S.-H. Kim, V. N. Manoharan, D. J. Pine, and S.-M. Yang, J. Am. Chem. Soc., 127, 15968–15975 (2005). Self-organization of bidisperse colloids in water droplets.
  5. F. W. Starr and F. Sciortino, J. Phys.: Condens. Matter, jul (2006). Model for assembly and gelation of four-armed dna dendrimers.
  6. Z. Nie, D. Fava, E. Kumacheva, S. Zou, G. C. Walker, and M. Rubinstein, Nat. Mater., 08 (2007). Self-assembly of metal-polymer analogues of amphiphilic triblock copolymers.
  7. J. P. K. Doye, A. A. Louis, I.-C. Lin, L. R. Allen, E. G. Noya, A. W. Wilber, H. C. Kok, and R. Lyus, Phys. Chem. Chem. Phys., 9, 2197–2205 (2007). Controlling crystallization and its absence: Proteins, colloids and patchy models.
  8. C. Mirkin, R. Letsinger, R. Mucic, and J. Storhoff., Nature, 382, 607 (1996). A dna-based method for rationally assembling nanoparticles into macroscopic materials.
  9. V. Workum and J. Douglas, Phys. Rev. E, 73, 031502 (2006). Symmetry, equivalence, and molecular self-assembly.
  10. A. Khan, Current Opinion in Colloid & Interface Science, 1, 614 (1996). Phase science of surfactants.
  11. P. van der Schoot and M. Cates, Langmuir, 10, 670–679 (1994). Growth, static light scattering, and spontaneous ordering of rodlike micelles.
  12. D. M. Kuntz and L. M. Walker, Soft Matter, 4, 286–293 (2008). Nematic phases observed in amphiphilic polyelectrolyte-surfactant aggregate solutions.
  13. J.-M. Jung and R. Mezzenga, Langmuir, 26(1), 504–514 (2010). Liquid crystalline phase behavior of protein fibers in water: Experiments versus theory.
  14. C. F. Lee, Phys. Rev. E, Sep (2009). Isotropic-nematic phase transition in amyloid fibrilization.
  15. A. Ciferri, Liq. Cryst., 34(6), 693–696 (2007). On collagen II fibrillogenesis.
  16. A. Aggeli, M. Bell, L. M. Carrick, C. W. G. Fishwick, R. Harding, P. J. Mawer, S. E. Radford, A. E. Strong, and N. Boden, J. Am. Chem. Soc., 125(32), 9619–9628 (2003). ph as a trigger of peptide -sheet self-assembly and reversible switching between nematic and isotropic phases.
  17. C. Robinson, Tetrahedron, 13(1-3), 219 – 234 (1961). Liquid-crystalline structures in polypeptide solutions.
  18. F. Livolant, A. M. Levelut, J. Doucet, and J. P. Benoit, Nature, 06 (1989). The highly concentrated liquid-crystalline phase of dna is columnar hexagonal.
  19. K. Merchant and R. L. Rill, Biophys. J., 12 (1997). Dna length and concentration dependencies of anisotropic phase transitions of dna solutions.
  20. F. Tombolato and A. Ferrarini, J. Chem. Phys., 122(5), 054908 (2005). From the double-stranded helix to the chiral nematic phase of b-dna: A molecular model.
  21. F. Tombolato, A. Ferrarini, and E. Grelet, Phys. Rev. Lett., Jun (2006). Chiral nematic phase of suspensions of rodlike viruses: Left-handed phase helicity from a right-handed molecular helix.
  22. E. Barry, D. Beller, and Z. Dogic, Soft Matter, 5, 2563–2570 (2009). A model liquid crystalline system based on rodlike viruses with variable chirality and persistence length.
  23. E. Grelet and S. Fraden, Phys. Rev. Lett., May (2003). What is the origin of chirality in the cholesteric phase of virus suspensions?
  24. S. Tomar, M. M. Green, and L. A. Day, J. Am. Chem. Soc., 129(11), 3367–3375 (2007). Dna-protein interactions as the source of large-length-scale chirality evident in the liquid crystal behavior of filamentous bacteriophages.
  25. A. Minsky, E. Shimoni, and D. Frenkiel-Krispin, Nat. Rev. Mol. Cell. Biol., 01 (2002). Stress, order and survival.
  26. J. Lydon, J. Mater. Chem., 20, 10071–10099 (2010). Chromonic review.
  27. K. Liu, Z. Nie, N. Zhao, W. Li, M. Rubinstein, and E. Kumacheva, Science, 329(5988), 197–200 (2010). Step-growth polymerization of inorganic nanoparticles.
  28. M. Nakata, G. Zanchetta, B. D. Chapman, C. D. Jones, J. O. Cross, R. Pindak, T. Bellini, and N. A. Clark, Science, 318, 1276– (2007). End-to-end stacking and liquid crystal condensation of 6 to 20 base pair dna duplexes.
  29. G. Zanchetta, M. Nakata, M. Buscaglia, N. A. Clark, and T. Bellini, J. Phys.: Condens. Matter, 20(49), 494214 (2008). Liquid crystal ordering of dna and rna oligomers with partially overlapping sequences.
  30. G. Zanchetta, F. Giavazzi, M. Nakata, M. Buscaglia, R. Cerbino, N. A. Clark, and T. Bellini, Proc. Natl. Acad. Sci. USA, 107(41), 17497–17502 (2010). Right-handed double-helix ultrashort dna yields chiral nematic phases with both right- and left-handed director twist.
  31. K. M. Guckian, B. A. Schweitzer, R. X.-F. Ren, C. J. Sheils, D. C. Tahmassebi, and E. T. Kool, J. Am. Chem. Soc., 122(10), 2213–2222 (2000). Factors contributing to aromatic stacking in water: evaluation in the context of dna.
  32. T. Bellini, R. Cerbino, and G. Zanchetta in Liquid Crystals - Materials Design and Self-Assembly, C. Tschierske, Ed., Vol. 318 of Topics in Current Chemistry; Springer Berlin / Heidelberg, 2012; pages 225–279.
  33. T. Kouriabova, M. Betterton, and M. Glaser, J. Mat. Chem., 20, 10366–10383 (2010). Linear aggregation and liquid-crystalline order: comparison of monte carlo simulation and analytic theory.
  34. X. Lü and J. Kindt, J. Chem. Phys., 120, 10328–10338 (2004). Monte carlo simulation of the self-assembly and phase behavior of semiflexible equilibrium polymers.
  35. G. J. Vroege and H. N. W. Lekkerkerker, Rep. Prog. Phys., 55, 1241–1309 (1992). Phase transitions in lyotropic colloidal and polymer liquid crystals.
  36. M. Dijkstra and D. Frenkel, Phys. Rev. E, Jun (1995). Simulation study of the isotropic-to-nematic transitions of semiflexible polymers.
  37. A. Khokhlov and A. Semenov, Physica, 108A, 546–556 (1981). Liquid-crystalline ordering in the solution of long persistent chains.
  38. A. Khokhlov and A. Semenov, Physica, 112A, 605–614 (1982). Liquid-crystalline ordering in the solution of partially flexible macromolecules.
  39. P. P. F. Wessels and B. M. Mulder, J. Phys.: Condens. Matter, 18(41), 9335 (2006). Isotropic-to-nematic transition in liquid-crystalline heteropolymers: I. formalism and main-chain liquid-crystalline polymers.
  40. M. Dennison, M. Dijkstra, and R. van Roij, Phys. Rev. Lett., May (2011). Phase diagram and effective shape of semiflexible colloidal rods and biopolymers.
  41. Z. Wang, D. Kuckling, and D. Johannsmann, Soft Mater., 1(3), 353–364 (2003). Temperature-induced swelling and de-swelling of thin poly(n-isopropylacrylamide) gels in water: Combined acoustic and optical measurements.
  42. Z. Y. Chen, Macromolecules, 26(13), 3419–3423 (1993). Nematic ordering in semiflexible polymer chains.
  43. T. Odijk, Macromolecules, 19, 2313–2329 (1986). Theory of lyotropic polymer liquid crystals.
  44. C. De Michele, T. Bellini, and F. Sciortino, Macromolecules, 45(2), 1090–1106 (2012). Self-assembly of bifunctional patchy particles with anisotropic shape into polymers chains: Theory, simulations, and experiments.
  45. C. De Michele, J. Comput. Phys., 229, 3276–3294 (2010). Simulating hard rigid bodies.
  46. F. Sciortino, C. De Michele, S. Corezzi, J. Russo, E. Zaccarelli, and P. Tartaglia, Soft Matter, 5, 2571—2575 (2009). A parameter-free description of the kinetics of formation of loop-less branched structures and gels.
  47. S. Corezzi, C. De Michele, E. Zaccarelli, P. Tartaglia, and F. Sciortino, J. Phys. Chem. B, 113(5), 1233–1236 (2009). Connecting irreversible to reversible aggregation: Time and temperature.
  48. M. Wertheim, J. Stat. Phys., 35, 19–34 (1984). Fluids with highly directional attractive forces. i. statistical thermodynamics.
  49. M. Wertheim, J. Stat. Phys., 35, 35–47 (1984). Fluids with highly directional attractive forces. ii. thermodynamic perturbation theory and integral equations.
  50. M. Wertheim, J. Stat. Phys., 42, 459–476 (1986). Fluids with highly directional attractive forces. iii. multiple attraction sites.
  51. L. Onsager, Ann. N. Y. Acad. Sci., page 627 (1949). The effects of shape on the interaction of colloidal particles.
  52. T. E. Ouldridge, A. A. Louis, and J. P. K. Doye, J. Chem. Phys., 134(8), 085101 (2011). Structural, mechanical, and thermodynamic properties of a coarse-grained dna model.
  53. T. E. Ouldridge, A. A. Louis, and J. P. K. Doye, Phys. Rev. Lett., 104, 178101 (2010). Dna nanotweezers studied with a coarse-grained model of dna.
  54. J. SantaLucia, Proc. Natl. Acad. Sci. USA, 95(4), 1460 (1998). A unified view of polymer, dumbbell, and oligonucleotide dna nearest-neighbor thermodynamics.
  55. J. A. Holbrook, M. W. Capp, R. M. Saecker, and M. T. Record, Biochemistry, 38(26), 8409–8422 (1999). Enthalpy and heat capacity changes for formation of an oligomeric dna duplex: Interpretation in terms of coupled processes of formation and association of single-stranded helices.
  56. T. E. Ouldridge Coarse-Grained Modelling of DNA and DNA self-assembly PhD thesis, University of Oxford, (2011).
  57. E. Protozanova, P. Yakovchuk, and M. Frank-Kamenetskii, J. Mol. Biol., 09 (2004). Stacked-unstacked equilibrium at the nick site of dna.
  58. P. Yakovchuk, E. Protozanova, and M. D. Frank-Kamenetskii, Nucl. Acids Res., 34(2), 564–574 (2006). Base-stacking and base-pairing contributions into thermal stability of the dna double helix.
  59. S. Bommarito, N. Peyret, and J. J. SantaLucia, Nucl. Acids Res., 28(9), 1929–1934 (2000). Thermodynamic parameters for dna sequences with dangling ends.
  60. L. Rovigatti and F. Sciortino, Mol. Phys., 109(23-24), 2889–2896 (2011). Self and collective correlation functions in a gel of tetrahedral patchy particles.
  61. M. Allen and G. Germano, Mol. Phys., 104(20-21), 3225–3235 (2006). Expressions for forces and torques in molecular simulations using rigid bodies.
  62. F. Romano, P. Šulc, L. Rovigatti, T. E. Ouldridge, A. A. Louis, and J. P. K. Doye (2012). In preparation.
  63. J. Parsons, Phys. Rev. A, 19, 1225–1230 (1979). Nematic ordering in a system of rods.
  64. C. Maffeo, B. Luan, and A. Aksimentiev, Nucl. Acids Res. (2012). End-to-end attraction of duplex dna.
  65. F. Sciortino, E. Bianchi, J. F. Douglas, and P. Tartaglia, J. Chem. Phys., 126, 194903 (2007). Self-assembly of patchy particles into polymer chains: A parameter-free comparison between wertheim theory and monte carlo simulation.
  66. E. Frezza, F. Tombolato, and A. Ferrarini, Soft Matter, 7, 9291–9296 (2011). Right- and left-handed liquid crystal assemblies of oligonucleotides: phase chirality as a reporter of a change in non-chiral interactions?
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description