Selfassembly of short DNA duplexes: from a coarsegrained model to experiments through a theoretical link
Abstract
Short bluntended DNA duplexes comprising 6 to 20 base pairs selfassemble into polydisperse semiflexible 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 selfassembly of such doublehelical duplexes with a combined numerical and theoretical approach. We simulate the bulk system employing the coarsegrained 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 parameterfree theoretical predictions provide an accurate description of the simulation results in the isotropic phase. In addition, the theoretical isotropicnematic phase boundaries are in line with experimental findings, providing a route to estimate the stacking free energy.
I Introduction
Selfassembly 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 selfassembly ubiquitous in nature and of interest in several fields, including material science, soft matter and biophysics (1); (2); (3); (4); (5). Through selfassembly 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 selfassembly 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 Bform 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 selfassembly 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 selfassembly of DNA duplexes provides a suitable way to access and quantify hydrophobic coaxial stacking interactions.
In order to extract quantitative informations from DNADNA coaxial stacking experiments, reliable computational models and theoretical frameworks are needed. Recent theoretical approaches have focused on the isotropicnematic (IN) transition in selfassembling systems (33); (34), building on previous work on rigid and semiflexible 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 superquadrics with quasicylindrical 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 coarsegrained model of DNA recently proposed by Ouldridge et al. (52), where nucleotides are modeled as rigid bodies interacting with sitesite 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 isotropicnematic phase boundaries.
To validate the theoretical predictions, we perform largescale 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 parameterfree 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 coarsegrained 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
We implement a coarsegrained model for DNA recently developed by Ouldridge et al. (53); (52). In such model, designed via a topdown 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 Bform. 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, WatsonCrick hydrogen bonding, stacking, crossstacking and coaxialstacking. 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 14base oligomers by Holbrook et al. (55). Hydrogenbonding and crossstacking potentials were adjusted to give duplex and hairpin formation thermodynamics consistent with the SantaLucia parameterization of the nearestneighbor 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 WatsonCrick 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 coaxialstacking interaction acts between any two nonbonded nucleotides and is responsible for the duplexduplex 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 CPUGPU code used in a previous work (60), and extend it to support the forcefields (61). Harvesting the power of modern Graphical Processing Units (GPUs) results in a fold speedup. The CPU version of the code, as well as the Python library written to simplify generation of initial configurations and postprocessing 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. cylinderlike 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 .
Iii Theory
We build on the theoretical framework previously developed to account for the linear aggregation and collective ordering of quasicylindrical 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
(1)  
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 ParsonsLee factor (63)
(2) 
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):
(3) 
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
(4) 
where
(5) 
With this choice for the free energy in Eq. (1) becomes:
(6)  
Minimization of the free energy in Eq. (6) with respect to provides the following expression for the average chain length :
(7) 
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.:
(8) 
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:
(9) 
where the term is the endmidsection contribution to the excluded volume of two hard cylinders (see Appendix B) and
(10) 
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:
Assuming that the orientational entropy can be approximated with the expression valid for long chains (43), minimization with respect to results in
(12) 
while using the approximated expression for short chains (43), one obtains
(13) 
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:
(14) 
Further refinements of the theory could be obtained by including a more accurate description of the orientational distribution in the proximity of the IN 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:
(15) 
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 largescale 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:
(16) 
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 endtoend bonded duplexes are not wellknown 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 largescale fullatom simulations by Maffeo et al. (64). They found that bluntended 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 .
In the present model the continuity of the helix under endtoend 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 BDNA 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 endtoend 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)
(17) 
where (65)
(18) 
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:

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

Set the counter and the energy factor .

Choose randomly two configurations and from the generated ensemble.

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 .

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,
(19) 
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:
(20) 
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
(21) 
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 .
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.
iv.4 Phase Coexistence: Theoretical predictions
A numerical evaluation of the phase coexistence between the isotropic and the nematic phases for the coarsegrained model adopted in this study is still impossible to obtain given the current computational power. We thus limit ourselves to the evaluation of the IN 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 IN 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 isotropicnematic coexistence region.
iv.5 Comparison between theory and experiments
The theoretical predictions concerning the isotropicnematic coexisting concentrations can be compared to the experimental results reported in Refs. (28) and (30) for bluntended 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 bluntend 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 endtoend bonds with more than one preferred azimuthal angle would increase the entropy of bonding, thus effectively lowering . Allowing for both left and righthanded 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 IN 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 backbonebackbone interaction of our coarsegrained 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.
iv.6 Comparison with Onsager Theory
The experimental average aggregation numbers are estimated in Refs. (13); (28) by mapping the selfassembled system onto an “equivalent” monodisperse system of hard rods with an aspect ratio equal to . In Ref. (44) it has been shown that the theoretically estimated isotropicnematic 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 reevaluated 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 reentrant behavior of the transition lines in the plane is observed. The reentrant behavior occurs for values of the stacking free energy accessed at temperatures between and . We believe that the reentrancy of the transition lines in the plane is a peculiar mark of the system polydispersity resulting from the reversible selfassembling into chains, and as such it should be also observable experimentally.
V Conclusions
In this article, we have provided the first study of a bulk solution of blunt ended DNA duplexes undergoing reversible selfassembly 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 selfassembly 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 IN 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 singlepeaked, is too restraining. In this respect, the present study call for an improvement of the coarsegrained 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 coarsegrained DNA model. Such procedure shows that values of the stacking free energy between and are compatible with the experimental location of the IN 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 allatom 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 allatom forcefield 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 monodisperse thin rods (28); (13). We have found that such approximation significantly underestimates at the IN transition concentration . In addition, the theoretical approach predicts a reentrant behavior of the transition lines in the  plane, a distinct feature of the polydisperse nature of the equilibrium chains.
Acknowledgments
We thank Thomas Ouldridge, Flavio Romano and Teun Vissers for fruitful discussions. We acknowledge support from ERC226207PATCHYCOLLOIDS and ITN234810COMPLOIDS 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 :
(22) 
where the functions , , describe the angular dependence of the excluded volume. The orientational probability is normalized such that
(23) 
The three contributions to the excluded volume in Eq. (22) come from endend, endmidsection and midsectionmidsection 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
(24) 
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.
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
(25)  
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 twodimensional fit to numerical data for using Eq. (9) as fitting function.
Appendix B Excluded volume of hard cylinders
For two rigid chains of length and which are composed of hard cylinders (HCs) of diameter and length , can be described by
(26)  
where , and are the orientations of two HCs and is the complete elliptical integral
(27) 
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)
(28) 
where
(29)  
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:
(30)  
(31)  
(32) 
The coefficient values resulting from the fitting procedure are , , , , , and .
Footnotes
 To whom correspondence should be addressed. Tel: +390649913524; Fax: +39064463158; Email: cristiano.demichele@roma1.infn.it
References
 I. Hamley, Introduction to Soft Matter, Wiley & Sons, 2007.
 S. C. Glotzer, Science, 306, 419–420 (2004). Some assembly required.
 G. M. Whitesides and M. Boncheva, Proc. Natl. Acad. Sci. USA, 99(8), 4769–4774 (2002). Beyond molecules: Selfassembly of mesoscopic and macroscopic components.
 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). Selforganization of bidisperse colloids in water droplets.
 F. W. Starr and F. Sciortino, J. Phys.: Condens. Matter, jul (2006). Model for assembly and gelation of fourarmed dna dendrimers.
 Z. Nie, D. Fava, E. Kumacheva, S. Zou, G. C. Walker, and M. Rubinstein, Nat. Mater., 08 (2007). Selfassembly of metalpolymer analogues of amphiphilic triblock copolymers.
 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.
 C. Mirkin, R. Letsinger, R. Mucic, and J. Storhoff., Nature, 382, 607 (1996). A dnabased method for rationally assembling nanoparticles into macroscopic materials.
 V. Workum and J. Douglas, Phys. Rev. E, 73, 031502 (2006). Symmetry, equivalence, and molecular selfassembly.
 A. Khan, Current Opinion in Colloid & Interface Science, 1, 614 (1996). Phase science of surfactants.
 P. van der Schoot and M. Cates, Langmuir, 10, 670–679 (1994). Growth, static light scattering, and spontaneous ordering of rodlike micelles.
 D. M. Kuntz and L. M. Walker, Soft Matter, 4, 286–293 (2008). Nematic phases observed in amphiphilic polyelectrolytesurfactant aggregate solutions.
 J.M. Jung and R. Mezzenga, Langmuir, 26(1), 504–514 (2010). Liquid crystalline phase behavior of protein fibers in water: Experiments versus theory.
 C. F. Lee, Phys. Rev. E, Sep (2009). Isotropicnematic phase transition in amyloid fibrilization.
 A. Ciferri, Liq. Cryst., 34(6), 693–696 (2007). On collagen II fibrillogenesis.
 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 selfassembly and reversible switching between nematic and isotropic phases.
 C. Robinson, Tetrahedron, 13(13), 219 – 234 (1961). Liquidcrystalline structures in polypeptide solutions.
 F. Livolant, A. M. Levelut, J. Doucet, and J. P. Benoit, Nature, 06 (1989). The highly concentrated liquidcrystalline phase of dna is columnar hexagonal.
 K. Merchant and R. L. Rill, Biophys. J., 12 (1997). Dna length and concentration dependencies of anisotropic phase transitions of dna solutions.
 F. Tombolato and A. Ferrarini, J. Chem. Phys., 122(5), 054908 (2005). From the doublestranded helix to the chiral nematic phase of bdna: A molecular model.
 F. Tombolato, A. Ferrarini, and E. Grelet, Phys. Rev. Lett., Jun (2006). Chiral nematic phase of suspensions of rodlike viruses: Lefthanded phase helicity from a righthanded molecular helix.
 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.
 E. Grelet and S. Fraden, Phys. Rev. Lett., May (2003). What is the origin of chirality in the cholesteric phase of virus suspensions?
 S. Tomar, M. M. Green, and L. A. Day, J. Am. Chem. Soc., 129(11), 3367–3375 (2007). Dnaprotein interactions as the source of largelengthscale chirality evident in the liquid crystal behavior of filamentous bacteriophages.
 A. Minsky, E. Shimoni, and D. FrenkielKrispin, Nat. Rev. Mol. Cell. Biol., 01 (2002). Stress, order and survival.
 J. Lydon, J. Mater. Chem., 20, 10071–10099 (2010). Chromonic review.
 K. Liu, Z. Nie, N. Zhao, W. Li, M. Rubinstein, and E. Kumacheva, Science, 329(5988), 197–200 (2010). Stepgrowth polymerization of inorganic nanoparticles.
 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). Endtoend stacking and liquid crystal condensation of 6 to 20 base pair dna duplexes.
 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.
 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). Righthanded doublehelix ultrashort dna yields chiral nematic phases with both right and lefthanded director twist.
 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.
 T. Bellini, R. Cerbino, and G. Zanchetta in Liquid Crystals  Materials Design and SelfAssembly, C. Tschierske, Ed., Vol. 318 of Topics in Current Chemistry; Springer Berlin / Heidelberg, 2012; pages 225–279.
 T. Kouriabova, M. Betterton, and M. Glaser, J. Mat. Chem., 20, 10366–10383 (2010). Linear aggregation and liquidcrystalline order: comparison of monte carlo simulation and analytic theory.
 X. Lü and J. Kindt, J. Chem. Phys., 120, 10328–10338 (2004). Monte carlo simulation of the selfassembly and phase behavior of semiflexible equilibrium polymers.
 G. J. Vroege and H. N. W. Lekkerkerker, Rep. Prog. Phys., 55, 1241–1309 (1992). Phase transitions in lyotropic colloidal and polymer liquid crystals.
 M. Dijkstra and D. Frenkel, Phys. Rev. E, Jun (1995). Simulation study of the isotropictonematic transitions of semiflexible polymers.
 A. Khokhlov and A. Semenov, Physica, 108A, 546–556 (1981). Liquidcrystalline ordering in the solution of long persistent chains.
 A. Khokhlov and A. Semenov, Physica, 112A, 605–614 (1982). Liquidcrystalline ordering in the solution of partially flexible macromolecules.
 P. P. F. Wessels and B. M. Mulder, J. Phys.: Condens. Matter, 18(41), 9335 (2006). Isotropictonematic transition in liquidcrystalline heteropolymers: I. formalism and mainchain liquidcrystalline polymers.
 M. Dennison, M. Dijkstra, and R. van Roij, Phys. Rev. Lett., May (2011). Phase diagram and effective shape of semiflexible colloidal rods and biopolymers.
 Z. Wang, D. Kuckling, and D. Johannsmann, Soft Mater., 1(3), 353–364 (2003). Temperatureinduced swelling and deswelling of thin poly(nisopropylacrylamide) gels in water: Combined acoustic and optical measurements.
 Z. Y. Chen, Macromolecules, 26(13), 3419–3423 (1993). Nematic ordering in semiflexible polymer chains.
 T. Odijk, Macromolecules, 19, 2313–2329 (1986). Theory of lyotropic polymer liquid crystals.
 C. De Michele, T. Bellini, and F. Sciortino, Macromolecules, 45(2), 1090–1106 (2012). Selfassembly of bifunctional patchy particles with anisotropic shape into polymers chains: Theory, simulations, and experiments.
 C. De Michele, J. Comput. Phys., 229, 3276–3294 (2010). Simulating hard rigid bodies.
 F. Sciortino, C. De Michele, S. Corezzi, J. Russo, E. Zaccarelli, and P. Tartaglia, Soft Matter, 5, 2571—2575 (2009). A parameterfree description of the kinetics of formation of loopless branched structures and gels.
 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.
 M. Wertheim, J. Stat. Phys., 35, 19–34 (1984). Fluids with highly directional attractive forces. i. statistical thermodynamics.
 M. Wertheim, J. Stat. Phys., 35, 35–47 (1984). Fluids with highly directional attractive forces. ii. thermodynamic perturbation theory and integral equations.
 M. Wertheim, J. Stat. Phys., 42, 459–476 (1986). Fluids with highly directional attractive forces. iii. multiple attraction sites.
 L. Onsager, Ann. N. Y. Acad. Sci., page 627 (1949). The effects of shape on the interaction of colloidal particles.
 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 coarsegrained dna model.
 T. E. Ouldridge, A. A. Louis, and J. P. K. Doye, Phys. Rev. Lett., 104, 178101 (2010). Dna nanotweezers studied with a coarsegrained model of dna.
 J. SantaLucia, Proc. Natl. Acad. Sci. USA, 95(4), 1460 (1998). A unified view of polymer, dumbbell, and oligonucleotide dna nearestneighbor thermodynamics.
 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 singlestranded helices.
 T. E. Ouldridge CoarseGrained Modelling of DNA and DNA selfassembly PhD thesis, University of Oxford, (2011).
 E. Protozanova, P. Yakovchuk, and M. FrankKamenetskii, J. Mol. Biol., 09 (2004). Stackedunstacked equilibrium at the nick site of dna.
 P. Yakovchuk, E. Protozanova, and M. D. FrankKamenetskii, Nucl. Acids Res., 34(2), 564–574 (2006). Basestacking and basepairing contributions into thermal stability of the dna double helix.
 S. Bommarito, N. Peyret, and J. J. SantaLucia, Nucl. Acids Res., 28(9), 1929–1934 (2000). Thermodynamic parameters for dna sequences with dangling ends.
 L. Rovigatti and F. Sciortino, Mol. Phys., 109(2324), 2889–2896 (2011). Self and collective correlation functions in a gel of tetrahedral patchy particles.
 M. Allen and G. Germano, Mol. Phys., 104(2021), 3225–3235 (2006). Expressions for forces and torques in molecular simulations using rigid bodies.
 F. Romano, P. Šulc, L. Rovigatti, T. E. Ouldridge, A. A. Louis, and J. P. K. Doye (2012). In preparation.
 J. Parsons, Phys. Rev. A, 19, 1225–1230 (1979). Nematic ordering in a system of rods.
 C. Maffeo, B. Luan, and A. Aksimentiev, Nucl. Acids Res. (2012). Endtoend attraction of duplex dna.
 F. Sciortino, E. Bianchi, J. F. Douglas, and P. Tartaglia, J. Chem. Phys., 126, 194903 (2007). Selfassembly of patchy particles into polymer chains: A parameterfree comparison between wertheim theory and monte carlo simulation.
 E. Frezza, F. Tombolato, and A. Ferrarini, Soft Matter, 7, 9291–9296 (2011). Right and lefthanded liquid crystal assemblies of oligonucleotides: phase chirality as a reporter of a change in nonchiral interactions?