text:
# References

Network of Hydrogen Bonds
as a Medium for DNA Interaction in Solvents

V.L. Golo, E.I. Kats††, and Yu.M. Yevdokimov†††

Department of Mechanics and Mathematics

Moscow University, Moscow 119899, Russia

†† Landau Institute for Theoretical Physics

Russian Academy of Sciences

Moscow, Russia

††† Institute of Molecular Biology

Russian Academy of Sciences

Moscow, Russia

Abstract

We suggest that the DNA molecules could form the cholesteric phase owing to an interaction mediated by the network of the hydrogen bonds (H-network) in the solvent. The model admits of the dependence of the optical activity of the solution on the concentration of the PEG, and the change in the sense of the cholesteric twist due to the intercalation by the daunomicyn. Using the experimental data for the cholesteric phase of the DNA dispersion, we obtain a rough estimate for the energy given by our model, and show that it should be taken into account as well as the energy due to the steric repulsion, van der Waals, and electrostatic forces, generally used for studying the DNA molecules. The elastic constant of the H-network generating the interaction between the DNA molecules is determined by the energy due to the proton’s vibration in the hydrogen bonds.

1. It is generally accepted that the formation of the cholesteric phases of the DNA dispersions is intimately related to the interplay among the basic properties of the DNA and the solvent, but so far the nature of the interaction, in the dispersions, between the molecules of DNA has remained obscure. In this communication we wish to suggest a model based on the structure of the network of the hydrogen bonds in water, which could explain the cholesteric transitions in the DNA dispersions. The main point is that the conformation of the network of the hydrogen bonds is influenced by the presence of the DNA molecules, so that the energy of the network is affected, and therefore one may expect a feed-back effect resulting in an interaction between the DNA molecules themselves.

It is worth noting that the hydrogen bonds, even though weak, $`5kcal/mol`$, that is much less than the covalent ones, are still larger than the van der Waals ones, and strong enough to hold the water molecules in a rigid structure in solid water (ice), . The structures of the water and the ice are similar, and either of them is determined by the network of the hydrogen bonds, , even though in liquid water some of the hydrogen bonds are broken; their number depending on temperature. The whole may be visualized as the fisherman’s net; the knots being atoms of the oxygen and the threads corresponding to the hydrogen bonds. The conformation of the net and the distribution of the protons are determined by Bernal-Fowler’s rules, i.e. there are two protons in the neighbourhood of an atom of the oxygen, and one proton is assigned to each hydrogen bond. We should allow also for defects; some of the threads (H-bonds) being ruptured owing to thermal fluctuations.

The conformation of the network of the hydrogen bonds in the DNA dispersions has certain properties that make it quite different from the network of the hydrogen bonds in the bulk solvent, so that the interaction between the DNA molecules at a distance less than $`50\AA `$ may result in new phenomena. The reason for this is twofold.

First, the quantum mechanical computations performed in a number of papers, , , show that the water molecules close to the DNA, i.e. at a distance of a few $`\AA `$, form a structure that follows the geometry of the molecule of DNA so that the set of the H-bonds is strictly determined by the DNA conformation, whereas in the bulk, at a distance of $`10\AA `$ to $`50\AA `$, the influence of the DNA is considerably weakend. Starting from classical crystallogrphic arguments, N.A. Buljenkov, , arrived at a conclusion that the main factor in the making of a conformation for the ensemble of the water and the DNA molecules, is the H-bonds. According to , inside one of the grooves of the DNA the water molecules form a kind of helix, whereas outiside the remaining water molecules are arranged in a hierarchical fashion, or, putting it in more quantitative terms, they should form a fractal structure according to a recursive algorithm, .

Second, the numerical simulation by Paola Gallo, , indicates that there are two distinct sets of water molecules confined in a silica-pore, i.e. those sufficiently close to the wall, and those outside in the bulk. The molecules close to the surface of the pore tend to form H-bonds with the atoms of the substrate so that the H-network close to the walls is strongly distorted, whereas the water in the bulk region, a few layers outside of the surface, shows the boson peak anomally, which is generally considered as a precursor of the glass transition, . The size of the silica pore is similar to that of the inter-molecular distance of the DNA dispersions, so that allowing for the difference in the potentials of surface interaction, one may suggest that the water molecules confined in the DNA dispersion may exhibit a similar behaviour. In this respect it is worthwhile to recall papers , , that claim the importance of what they call hydration forces, which dominate at close separations between the DNA molecules, so that the resultant interaction resemble neither the pure electrostatic nor the steric one.

There are experimental data suggesting that there be a need for taking into account the part played by the inter-molecular medium of the DNA dispersions. Indeed, the transitions from the cholesteric to the nematic phases, which are observed in the DNA dispersions , have a number of features that run considerable difficulties in explaining within the conventional framework of electrostatic and van der Waals forces, (see paper  in which the problem is studied within the framework of the electrostatic interaction between the molecules of the DNA).

Let us recall the most prominent experimental facts.

(i) The liquid crystalline phases formed from the double-stranded DNA, exhibit a negative band in the CD spectra, whereas the liquid crystalline-phases and dispersions formed from the RNA, under the same conditions, have a positive band in their CD spectra, . The positive (resp. negative) sense of the CD band is alleged to correspond to the left (resp. right) twist of cholesteric. Under the experimental conditions of , the difference between the right-helix of the B-form of the DNA and the right double-helix form of the RNA appears to be not sufficiently significant, there is only the change of thymine for uracil, to cause a substantial spatial effect solely because of the electrostatical and the van der Waals forces. In fact, one should expect only a tiny modification of the inter-molecular interaction, and not the reconstruction of the total structure.

(ii) In paper  it is reported that nucleic acid molecules form the cholesteric phases in the presence of the poly(ethylene glycol), (PEG), only at a certain osmotic pressure of the PEG. The osmotic pressure can be controlled through the variation of the concentration of PEG in the solvent; to a certain extent, it is similar to the classical experiment with a semipermeable piston that exerts a pressure equal to the osmotic pressure of the PEG solution. It is generally considered to be an efficient means for estimating the repulsion interaction between the DNA helices in solution, , . The results of paper  are illustrated in Fig.1; it appears that there should be a change in the state of the dispersion. There is an obvious correlation between the character of the osmotic curves and the optical activity of the solution, which should correspond to the formation of the cholesteric phase. But the osmotic effect is due to the action of the hydration force associated with the removal of the water molecules from the DNA environment. Therefore, the force generates the principal interaction between DNA molecules, rather than the van der Waals or the electrostatic one (cf. ).

(iii) The CD spectra of the liquid crystalline dispersions formed from the right handed-helices of the B-DNA has the important characteristic, the intensive negative band, whereas for the case of the A-DNA, the corresponding band is positive. One may change the CD-spectrum of the DNA dispersion by dyes, so that two bands instead of one appear. Depending on the dye one employs, the two bands may vary. If the dye is intercalated between the DNA base-pairs, the two bands are related to (1) the absorption region of the DNA nitrogen bases, and (2) the absorption region of the compound itself, respectively. If, as is the case with the daunomycin, a part of the compound lies ”outside”, i.e. in one of the grooves of the DNA, there appear two positive bands in the CD spectrum of the molecule. It is alleged that the two positive bands are the signature of the right-handed twist of the cholesteric helix formed by the DNA molecules. Thus, the sense of the twist is changed by the daunomycin. But the daunomycin can be eluted by sodium dodecylsulfate (SDS), so that the positive band related to the daunomycin disappears; nonetheless the positive band in the absorption band of the DNA still remains. One may suggest that the right-handed twist of the cholesteric structure of the DNA dispersion remains as well. Thus, it is possible to obtain the right-handed twist cholesteric structure, besides the usual left one, for the same B-form of the DNA. One may suggest that if the cholesteric twist were determined by the van der Waals and the electrostatic forces, as well as the form of the DNA, the co-existence of these two metastable states should be unlikely.

(iv) According to , even small local defects in the structure of the DNA molecules may cause the change in the cholesterical structure. It is important that defects of the surface of the DNA molecule result in the formation of a helical structure in which the sense of the cholesteric twist differs from the sense of the twist of the DNA cholesteric, and the intense band in the CD spectrum must tend to zero as the number of the defects increases. It is possible that the defects affect the neighborhood of a DNA molecule, thus preventing the formation of the cholesteric packing; the steric, van der Waals, and electrostatic effects not being involved, .

2. As was explained above, in the presence of the molecules of DNA, e.g. in the dispersions, the network of the hydrogen bonds has different shapes in the regions close to the DNA molecules and in those of the bulk, a few layers of the water molecules away from them. One may imagine that close to the DNA molecules, the net looks like being rolled round a pole; it is spreading more leisurely away from the DNA molecules.

The key problem is to find an estimate for the elastic properties of the system. In fact, the differences in energy between microscopical states of the network of the hydrogen bonds are small, and therefore the configuration of the system is determined by entropy. Therefore, the H-network may be considered, from the point of view of its elastic properties, as a system of the ”athermal” type, . Presently, the information about their elastic properties is scarce, so that one may expect that the study of the cholesteric phases of the DNA dispersions may provide a new insight into the matter (see below). We shall assume that the network of the H-bonds in solvent has a certain amount of elasticity, and in fact enough to account for the interaction of the DNA molecules. At the same time the system may also have some elastic energy of the usual thermal type; presumably, not much, depending on the elasticity of the H-bonds, i.e. their being extensible. The total energy comprized of the elastic and the entropy term, under constraints imposed by the boundary effects, should result in an effective iteraction between the DNA molecules, through the medium of the network of the hydrogen bonds. It is to be noted that we have neglected the van der Waals and the electrostatic forces, usually suggested to be of primary importance, on the ground of the arguments given in the preceeding section. Their effect could still be important in the system’s choosing its orientation as regards the cholesteric twist (see below), even though the determining factor should be the interaction mediated through the H-bonds.

The picture given above implies that there is a neighbourhood of the DNA molecule in the solvent in which the conformation of the set of the molecules of water is strongly influenced by the geometry of the DNA helix. This assumption is in full agreement with the theoretical and experimental results reported in the preceeding section. In this neighbourhood, which spreads out at a distance of approximately $`10\AA `$, the conformation of the H-bonds follows the grooves of the DNA and has the shape resembling that of the helix. In this respect it differs considerably from that in the bulk, (compare the previous section, ,  ).

According to the numerical simulation of  the water molecules in the major groove form a regular framework comprized of hydrogen bonded pentagons that lie in the same plane as the base-pairs; in contrast, the water molecules in the minor groove form only a monolayer. Consequently, the conformation of the network of the H-bonds close to a molecule of the DNA could be compared, at least qualitativeley, to a helical, or winding, dislocation. The latter has the specific feature of being a kind of singularity in the fabric of the bulk network of the H-bonds, and in this sense one may consider it as a topological dislocation.

If we assume that the elastic, perhaps totally ”athermal”, or entropy generated, interaction between the dislocations of the network of the H-bonds is at work, we may suggest as well, within the usual framework of elasticity theory, that two parallel dislocations of the type interact according to the usual logarithmic law, i.e. the energy per unit length of the interaction of two parallel dislocations is given by the equation

$$E_{disl}=\frac{ϵe_1e_2}{2}\mathrm{ln}r$$
(1)
in which $`ϵ`$ is the elastic modulii, $`e_1,e_2`$ are the ”charges” determined by the winding of the net at the cores of the dislocations, that is the molecules of the DNA, and $`r`$ is the distance between the dislocations. The interaction is the repulsive one. According to the arguments given above, Eq.(1) gives the energy of the interaction of the two parallel molecules of the DNA, through the medium of the network of the H-bonds.

Next, let us consider two molecules of the DNA that are not parallel. Since the cholesteric phases in DNA dispersions are studied for segments of DNA of sufficiently low molecular weight, we may suppose that the molecules are of a length $`l`$, smaller than the persistent one, and we may visualize them as segments of two skew lines. For the sake of simplicity, and allowing only for a rough estimate of the energy of their interaction, we suppose as well that their centres be on the same line, at a distance $`d`$; the segments being at an angle $`\varphi `$ to each other. We shall use the planar approximation for the interaction, that is we assume that only bits of the molecules lying in the same planes do interact in accord with Eq.(1), so that the corresponding energy reads

$$F_{bits}=\frac{ϵ}{2}\mathrm{ln}\left[d^2+4\left(\frac{x}{l}\right)^2\mathrm{sin}^2(\varphi /2)\right]+ϵ\mathrm{ln}a$$
(2)
where $`x`$ is the position of the bits on the line of the molecule, and $`a`$ is the radius of the core of the dislocation, i.e. the radius of the DNA molecule. To derive the energy of the total interaction of a pair of the DNA molecules, we shall take the integral of the function defined by Eq.(2)

$$U_{disl}=_{l/2}^{+l/2}F_{bits}𝑑x$$
(3)
The energy given by Eq.(3) is only a part of the total expression for the interaction energy of the DNA molecules in solvent. It is necessary to allow for the isotropic pair interation, $`U_0`$, the conventional energy due to the van der Waals, electrostatic and steric forces, $`U_1`$, the anisotropic contribution related to the chiral state, $`U_2`$, so that the total expression has the form

$$W_\varphi =U_0+U_1+U_2+U_{disl}$$

It is important that the term $`U_0`$ does not depend on the angle $`\varphi `$ describing the mutual orientation of the molecules. The twist angle, $`\varphi `$, is small. Therefore, the term $`U_1`$ is of the second order in $`\varphi `$, and $`U_2`$ of the third,

$$U_1=A\varphi ^2,U_2=B\varphi ^3$$
(4)
In what follows we shall assume that the contributions due to $`U_1`$, which takes into account the combined influence of the steric, van der Waals, and electrostatic forces, is comparable with $`U_{disl}`$, whereas $`U_2`$ is a small correction. It is also necessary to take into account the energy due to the proton transfer, but within the approximation used in this paper it does not depend on the cholesteric twist, $`\varphi `$, and we may postpone considering the effect until Section 3.

On expanding $`U_{disl}`$ to within the fourth order in $`\varphi `$, and neglecting for a while the contribution due to $`U_2`$, we obtain the following expression for $`W_\varphi `$

$$W_\varphi =\lambda \mathrm{ln}\left(\frac{d}{a}\right)+\left[A\frac{\lambda }{3}\left(\frac{l}{d}\right)^2\right]\varphi ^2+\frac{\lambda }{2}\left(\frac{l}{d}\right)^4\varphi ^4$$

in which $`\lambda =lϵ`$. It is worth noting that for the experimental data, i.e. $`d`$ within $`30\AA `$ to $`50\AA `$, $`l`$ of the size $`500\AA `$, and $`\varphi `$ of the order $`0.01rad`$, we have the small parameter $`\kappa =l\varphi /d`$, and Eq.(S0.Ex2) is an expansion of $`W_\varphi `$ in $`\kappa `$. The minimization of $`W_\varphi `$ gives the value of the cholesteric twist angle

$$\varphi _{}^2=\frac{d^2}{l^4}\left(\frac{A}{\lambda }d^2\frac{l^2}{3}\right)$$
(5)
At first sight, the cholesteric phase appears to exist for all intermolecular distances $`dd_0=l\sqrt{(3A/\lambda )}`$. In fact, the situation is different. For one thing, for small values of $`d`$, say less $`20\AA `$, the asumption that the H-network be soft enough, breaks down, for the hydrogen bonds form a rigid structure due to the surface effects of the DNA molecules, and in this region of $`d`$ the nematic phase is observed, , , . For another, the value of $`d_0`$ is determined by the relative sizes of the elastic constants $`A`$ and $`\lambda `$. According to , , , the cholesteric phase is observed for the values $`d30\AA `$ and $`l500\AA `$, therefore, we may conclude that $`A`$ could be less than $`\lambda `$ at least by an order of magnitude, so that $`d_0`$ may be of the same order as $`l`$. Thus, Eq.(5) describes the cholesteric twist $`\varphi _{}`$, within the range of $`30\AA `$ to $`50\AA `$ for $`d`$. Our discussion is very crude, and qualitative, so that we may only infer from Eq.(5) that the twist angle, $`\varphi _{}`$, increases with $`d`$ in the region indicated above, and in this respect our result agrees with the formular for $`\varphi _{}`$

$$\varphi _{}\frac{\sqrt{d}}{l}$$

obtained in paper  within the framework of the electrostatic approach. But the region of $`d`$ is too small to make any asymptotic conclusions, especially that a significant dependence of $`\varphi _{}`$ on the spatial separation $`d`$ has not been observed, so far. It should be noted, as well, that paper  neglects the conclusions of Parsegian et al, , , that the hydration forces be of primary importance, and does not consider the PEG experiments.

It is important that below the threshold we have the two values for the twist angle, $`\varphi _+`$ and $`\varphi _{}`$, owing to the energy $`U_2`$, (see Eq.(4)), which takes off the degeneracy in $`\pm \varphi _{}`$. Consequently, we may expect the existence of two metastable cholesteric structures, above the nematic threshold. This circumstance may have a bearing upon the change of the CD spectrum due to the daunomycin, described above. We suggest that under the action of the latter, in a layer of thickness of approximately $`5\AA `$, the conformation of the surface is deformed so that the structure of the hydrogen bonds in the whole volume is changed in such a way that the system comes down to the minimum different from that it would occupy in the absence of the daunomycin.

Conformational changes in the DNA molecules could be a cause for the reconstruction of the system of the hydrogen bonds in the bulk of the solvent, and therefore result in the formation of various cholesteric structures. In this respect, small local defects mentioned above, as well as the conformation structures of the double-helix of the B-DNA and the double helix form of the A-RNA, may bring about changes in the network of the hydrogen bonds that could result in different liquid crystalline structures. Similarly, the PEG solutions may provide conditions for the formation of different networks of the hydrogen bonds that could generate effective interactions between the DNA molecules that may result in liquid crystalline structures, The phenomena are difficult to explain within the framework of the theory using only electrostatic and van der Waals forces.

So far we have considered the network of the H-bonds generally with respect to its conformation, that is the knots corresponding to the positions of the atoms of oxygen and the threads to the H-bonds. It is important that stringent topological conditions should be imposed on the network owing to the requirements due to the structure of ice, , i.e. among $`16`$ possible combinations of the protons round an atom of the oxygen only $`6`$ are allowed, so that the residual entropy of ice turns out to be equal to $`S_0=Nk\mathrm{ln}(3/2)`$. The topological constraints are also essential for the proton dynamics of the H-bonds. In this respect it is worth noting the phenomenon reported in paper  that the protons could perform cooperative or concerted motions (flip-flop) by jumping from one state to another. Saenger et al, , claim that one could find chains of the water hydroxyl groups $`OH\mathrm{}OH\mathrm{}OH`$ in which the protons would oscillate clockwise and anti-clockwise in a flip-flop manner. Obviously, the choice of the circles allowing for the flip-flop motion is determined by the topology of the network of the H-bonds.

3. It is instructive to make a rough estimate of the effect produced by the proton’s vibration. One can visualize it as a motion parallel to the valence bond $`OH\mathrm{}`$, in a cell of the order of $`\sigma `$; $`\sigma `$ being the size of the hydrogen bond. In fact, the proton moves in the field of the double-well potential inside the hydrogen bond, with a characterisitic frequency, $`\nu _s`$, which is of the order $`3000cm^1`$; the corresponding energy being $`\mathrm{}\omega _s`$ per H-bond, . Here we take into account only the lowest frequency, and neglect the obertones. The number of these cells per area of the cross section of the volume in the DNA dispersion depends on the separation $`d`$ of the molecules. The value of $`\sigma `$ could not change appreciably with $`d`$. Consequently, the number of the hydrogen bonds per pair of the DNA molecules involved in the conformation of the cholesteric structure should be proportional to

$$N_H=\frac{l}{\sigma }\left(\frac{d}{\sigma }\right)^2$$
(6)
so that the energy due to the proton vibration may be estimated at

$$U_P=\frac{C}{2}\frac{l}{\sigma }\left(\frac{d}{\sigma }\right)^2\mathrm{}\omega _s$$
(7)
where $`C`$ is a constant, presumably of order $`1`$. It is obviously similar to the zero point energy, familiar in the theory of the so-called Casimir effect.

The zero point energy may play a significant role in determining macroscopical properties of a system. For example, its large value for the helium atoms is the reason why the volume of liquid helium at zero pressure is almost three times as large as the volume it would have had in its absence, . In the situation of the network of the hydrogen bonds, it is important that the energy of the proton’s vibration, $`\mathrm{}\omega _s\mathrm{6\hspace{0.17em}10}^{13}erg`$, is approximately $`16\%`$ of the energy of the hydrogen bond, $`5kcal/mol`$, or $`\mathrm{3.46\hspace{0.17em}10}^{12}erg`$. Hence, one may expect that the interplay among the elastic energy of the H-network and the vibrational energy of the protons should result in a meaningful structure. In fact, it does, and the expressions for $`W_\varphi `$ and $`U_P`$ having opposite signs insures that the formation of the dislocations in the H-network allows the system to lower its total energy. We may cast these arguments in a more quantitative form in the following way.

Using the expression for the energy $`W_\varphi `$ we may cast the energy of the interaction of the two molecules in the form

$$E_H=W_\varphi +U_P$$

and on minimizing $`E_H`$ we obtain Eq.(5) and the second equation, $`E/d=0`$, which can be cast in the form

$$\lambda C\frac{l}{\sigma }\left(\frac{d}{\sigma }\right)^2\mathrm{}\omega _s$$
(8)
From Eqs.(5-8), we may derive an estimate for the elastic constant $`\lambda `$. In fact, since $`\nu _s3000cm^1`$, and $`d`$ and $`l`$ of $`50\AA `$ and $`500\AA `$, respectively, (see ), we get $`\lambda 10^8erg`$. According to the arguments given above, $`A`$ is less than $`\lambda `$, by at least an order of magnitude.

At this point it is worthwhile to draw attention to the part played by the PEG. The latter is alleged to influence the hydration forces, i.e. to diminish the number of the H-bonds involved in the interaction of the DNA molecules. Therefore, the presence of the PEG should reduce the value of $`N_H`$ given by Eq.(6), and the constant $`C`$ should decrease as well. Consequently, in accord with Eq.(8), the elastic constant $`\lambda `$ should decrease and therefore the twist angle $`\varphi _{}`$ should diminish. But, the decrease in $`\varphi _{}`$ has not been observed. This is, perhaps, due to the circumstance that $`A`$ may significantly depend on $`d`$.

4. Our hypothesis as to the nature of the interaction of the DNA molecules that leads to the formation of the cholesteric phase, is in qualitative agreement with exprimental data; it accomodates

* the strong dependence of the formation of the cholesteric structure and their properties, on the conformational properties of the nucleic acids;
* the predominance of the hydration forces demonstrated with the experiments using the PEG solutions;
* the change in the sense of the cholesteric twist due to the intercalation by the daunomycin;
* the large effects caused by small defects of the molecular structure, and the different senses of the twist angles for the DNA and the RNA-molecules.

We have not found the isotope effect, and in fact so far it has not been observed.

But it is worth noting that the concept of the network of the H-bonds as a medium for the interaction between molecules of nucleic acids in dispersions, makes for studying the water itself. Indeed, according to our model, one could employ the DNA dispersions as a valuable probe into the structure of water, and in this respect it is very promising that we have obtained the reasonable elastic characteristics of the network of the hydrogen bonds.

The authors are thankful to N.A. Buljenkov, V.I. Ivanov, and G.G. Malenkov, for the useful discussions.

Figure Captions

Fig.1
The dependence of $`\mathrm{log}\mathrm{\Pi }`$ on average distance between nucleic acid molecular axes $`(D,\AA )`$.
(1) DNA;
(2) RNA.

Fig.2

A.
The CD spectra of the liquid-crystalline dispersion (curve 1) treated by DAU (curves 2,3).
$`C_{DNA}`$$`28.2\mu g/ml`$ ;
$`C_{PEG}`$$`150\mu g/ml`$ ;
$`(0.3MNaCl+0.005Naphosphatebuffer)`$ ;
Curve 1 — $`C_{DAU}=0`$;
Curve 2 — $`C_{DAU}=6.810^6M`$;
Curve 3 — $`C_{DAU}=27.210^6M`$;
$`\mathrm{\Delta }A`$ in $`mm`$; $`1mm`$$`10^5`$ opt. units.

B.
The CD spectra of the liquid-crystalline dispersions formed by (DNA-DAU) complexes.
$`C_{DNA}`$$`28.2\mu g/ml`$ ;
$`C_{PEG}`$$`150mg/ml`$ ;
$`(0.3MNaCl+0.005Naphosphatebuffer)`$ ;
Curve 1 — $`C_{DAU}=13.610^6M`$ ;
Curve 2 — $`C_{DAU}=37.610^6M`$ ;
$`\mathrm{\Delta }A`$ in $`mm`$; $`1mm`$$`10^5`$ opt. units.


format: text

text:
# The Dynamical Interstellar Medium: Insights from Numerical Models

## 1. Introduction

The modern picture of the interstellar medium (ISM) as a pervasive medium containing gas at temperatures ranging from hundreds to millions of degrees was first hinted at by Copernicus observations of the far ultraviolet absorption lines of O vi (Jenkins & Meloy 1974, Jenkins 1978). The ambient interstellar radiation field is too soft to photoionize O vi, so it must be thermally ionized by $`10^5`$ K gas. However, the cooling curve for gas with solar abundances peaks around $`10^5`$ K, implying that this gas exists in thin, cooling layers at the surface of a reservoir of $`10^6`$ K gas. This led to models of a thermally regulated, three-phase ISM controlled by the balance between energy input from supernovae and energy dissipation by radiative cooling (Cox & Smith 1974; McKee & Ostriker 1977). A new model is now beginning to emerge, in which dynamical compression rather than thermal conduction determines the thermal state of the gas as well as its structure.

The major question driving this new model is, what determines the star formation rate? The old idea that it might be determined by a balance between supersonic turbulence and self-gravity has come back into favor recently as the idea of magnetostatic support regulated by ambipolar diffusion has run into significant problems: magnetic fields weaker than required (Crutcher 1999); observed cores denser than predicted (Nakano 1998); and the inability of magnetic fields to maintain the observed supersonic motions (Mac Low et al. 1998, Stone, Ostriker, & Gammie 1998). In this paper I argue that numerical models point to supernovae as the main driver of that supersonic turbulence, at least in star-forming regions of the ISM.

As has long been recognized, the odd thing about the observed star formation rate in normal galaxies is how low it is. Typical molecular cloud free-fall times are of order 1 Myr, but there is still plentiful gas available after $`10^4`$ times that long. Another way to see this is to note that if the $`10^{10}M_{}`$ of gas observed in our galaxy were to collapse on timescales approaching the free-fall time, the star formation rate would be far higher than the observed $`12M_{}\text{ yr}^1`$.

This has led to the suggestion that molecular clouds are long-lived objects, with lifetimes ranging from 30 Myr (Blitz & Shu 1980) to as long as 100 Myr (Solomon & Sanders 1979, Scoville & Hersh 1979). In addition to the argument from the global star formation rate given above, Blitz & Shu (1980) also point to the extent of the presence of molecular clouds behind spiral arm shocks and the stellar ages then thought to be associated with the clouds as evidence for their lifetime estimate. Recently, however, Ballesteros-Paredes, Hartmann, & Vázquez-Semadeni (1999), have proposed much shorter lifetimes of only 5–7 Myr, based on the lack of post-T Tauri stars with ages in the 5–10 Myr range that are clearly associated with molecular clouds. Individual clouds might be quite short-lived, though cloud-forming regions might last significantly longer, explaining the observed spatial relationships with spiral arms.

## 2. Simulations

### 2.1. Questions

Large-scale simulations of the ISM address a number of questions in ways complementary to other methods. First, they can begin to reveal the topology of the different phases of the ISM, distinguishing between scenarios in which the hot gas is confined in isolated bubbles, is distributed in an interconnected tunnel network, or forms a uniform sea surrounding isolated clouds of warm and cold gas. The distribution of cold, neutral gas determines whether ionizing radiation from OB associations in the galactic plane can fully account for the ionization of the Reynolds layer far above it.

Second, the production of star-forming regions appears to be largely determined by large-scale flows in the ISM, influenced by the gravitational potential. Large-scale simulations can thus yield information on the formation and lifetimes of molecular clouds. Related issues can also be addressed, for example, whether supernovae and superbubbles are more important in triggering nearby star formation or in driving turbulence that supports the ISM against gravitational collapse and star formation.

Finally, two other major issues that I will not focus on to the same extent are the production of galactic magnetic fields in a turbulent dynamo, and the vertical flow of the ISM in galactic fountains and winds.

### 2.2. Objects vs. Interactions

Traditionally, numerical simulations of the ISM focussed on specific objects, such as supernova remnants, stellar wind bubbles, superbubbles, shocked clouds, or bow shocks produced by fast-moving stars or stellar jets. Such simulations have proved useful, but are limited by the lack of interactions with surrounding structures.

Observations of the ISM are structured at all scales. A frequent problem encountered by observers is distinguishing objects such as superbubbles from the structured background. The background structure is often as important in defining the final nature of an object as more traditional considerations such as the central energy source or cooling physics. The interactions between “objects” such as supernova remnants or superbubbles becomes as important as the objects themselves. The random supersonic flows generated by these interactions in turn determines whether regions will be supported against gravitational collapse, or be compressed and cooled, starting a process of collapse leading to molecular cloud and star formation.

### 2.3. Models

In the last five years, a significant number of global models have been computed. They are distinguished by the different physical processes included, and by the different numerical techniques used.

The most important distinction appears to be the driving mechanism for the modeled interstellar turbulence. The most abstract model for the driving is simply a uniform, Gaussian field, which has both the advantage and disadvantage of being very general, but not reproducing the specifics of particular real drivers. Three real drivers have also been used in various models: stellar ionization heating, rotational shear, and supernovae. The last appears to be the most realistic, although the other two add important physics when used in addition to the supernovae.

The dimensionality and coverage of the simulation are also very important. Two-dimensional simulations overstate the difficulty of hot gas escaping from the plane and cannot be used fruitfully to explore questions of magnetic field generation, but have yielded fundamental insights into the turbulent ISM. Three-dimensional models have been done in both local and global modes, either focussing on a particular region of a disk (typically around 1 kpc2) or attempting to encompass the entire disk.

Three main technical issues differentiate the models. First is the numerical diffusivity of the method used. Smoothed particle hydrodynamics (SPH) tends to be the most diffusive method due to the need to average over a fairly large number of particles to derive the flow quantities. Spectral methods are excellent for smooth flows, but handle shocks and other abrupt features in flows poorly. Grid-based methods have a numerical diffusivity determined by a combination of the order of the method, the advection technique used (e.g. second-order Van Leer, piecewise parabolic method (PPM), or total variation diminishing (TVD) methods), and the technique used for handling shocks.

The second main issue is the shock handling technique itself. Artificial viscosity can be used with particle, spectral, and grid methods to resolve shocks over several resolution lengths. Exact or linearized Riemann solvers resolve shocks more accurately and in fewer zones. They compute the flow at the edge of each zone by solving for the combination of shocks and compression or rarefaction waves that results from the breakup of the discontinuity at each zone boundary.

The third issue is the structure of the grid in those methods using one. The simplest grids are fixed, Eulerian grids, with periodic boundary conditions in the case of local simulations. Ratioed fixed grids that focus resolution toward a region of interest such as the galactic plane or the galactic center are the next step. SPH models effectively use an unstructured Lagrangian grid that follows the flow, with the advantage of high resolution in dense regions, but the corresponding disadvantage of low resolution in more rarefied regions. Adaptive mesh refinement focusses resolution into regions of interest regardless of density, at the cost of significant additional programming complexity.

## 3. Results

The outlines of a coherent picture of the interstellar medium is emerging from the numerical models, one that places a greater emphasis on dynamical interactions and less on equilibrium structures than previous analytically constructed pictures.

In the absence of supernova heating and the presence of normal interstellar cooling, the gas cools rapidly, forming dense clouds scattered through the disk of the galaxy that interact with each other, surrounded by a thicker disk of warm $`10^4`$ K gas, but without hot, $`10^6`$ K, X-ray emitting gas. An example of this solution is seen in the global SPH computations by Gerritsen & Icke (1997).

The global, two-dimensional PPM computations of Wada & Norman (1999) add the element of driving from galactic shear. Starting with a Toomre stable disk, they find that the gas rapidly cools and collapses into a filamentary network of clouds that are cold enough to be Toomre unstable. Although occasional cloud collisions generate small amounts of $`10^5`$ K gas, the interesting conclusion can be drawn from these computations that the Toomre stability criterion must be applied to a medium whose velocity dispersion is maintained by some mechanism other than ionization heating.

At the same time, it has become clear that dense interstellar clouds can be formed by turbulent compression. The local, 2D, MHD, spectral models of Passot, Vázquez-Semadeni, & Pouquet (1995) show that this can be a dominant mechanism for structure formation in the warm and cold ISM even in the presence of magnetic fields. Hennebelle & Pérault (1999, 2000) have explored the interplay of compression and thermal condensation due to conduction in detail with 1D, adaptive mesh models with and without magnetic fields, showing that compression will cause collapse on far shorter timescales than thermal condensation under typical ISM conditions. As thermal condensation is the dominant mechanism for cold cloud formation in analytic equilibrium models such as those of McKee & Ostriker (1977), this is an important conclusion.

Using supernovae as the main driving mechanism produces hot gas with moderate filling factor in the disk of galaxies with parameters similar to those of the Milky Way, and much larger filling factor above the plane. Local, 2D models done with a second-order, Van Leer algorithm on a ratioed grid by Rosen & Bregman (1995) established that such driving alone could lead to the formation of a network of cold clouds and hot gas, along with a flow of hot gas out of the disk. Use of adaptive mesh refinement with a Riemann solver allowed Avillez (1999) to perform 3D models extending 10 kpc into the upper halo. He finds that hot gas is not trapped in the disk as occurred in the 2D models. Ballesteros-Paredes, Avillez, & Mac Low (2000) show that instead, isolated, cold, dense clouds with masses and sizes typical of giant molecular clouds form in and close to the disk. Column density distributions of cool and cold gas from these models strongly resemble observations of H i, for example by Kim, Stavely-Smith & Bessell (1999).

These models either lack self-gravity, or, in the case of Passot et al. (1995) lack the resolution for self-gravity to become significant. However, isothermal, uniformly-driven models provide some insight into whether supernova-driven turbulence could support the gas against gravitational collapse. Klessen, Heitsch, & Mac Low (2000) show that supersonic turbulence can globally support a thermally Jeans unstable region, but local collapse will still occur due to the density enhancements caused by shock waves in the supersonic turbulence, except in the case of extremely strong turbulence. In fact, it appears that a signpost of global turbulent support under molecular cloud conditions is isolated star formation, while lack of turbulent support leads to efficient runaway star formation (Klessen, Burkert, & Bate 1998, Klessen & Burkert 2000a, 2000b). Heitsch, Mac Low, & Klessen (2000) show that magnetic fields add small-scale structure, and slow local collapse, but do not prevent it.

Applying the insights gained from these models to the interstellar case, it appears that the supernova-driven turbulence may naturally create regions that are gravitationally-bound, but then may drive enough turbulence within them to reduce their star-formation rate to the low levels observed in typical molecular clouds, probably with the help of interstellar magnetic fields. Observations of molecular clouds appear to show that they have self-similar structure to scales of as much as tens of pc (e.g. Stutzki et al. 1998), which Mac Low & Ossenkopf (2000) find to be most readily explained by driving of the turbulence from even larger scales, again consistent with large-scale driving by supernovae.

Korpi et al. (1999) use a 3D, fixed grid, local, MHD code to model the evolution of the magnetic field in shearing, supernova-driven turbulence. Korpi (1999) measures the stress terms in the flow to find preliminary evidence for the operation of a dynamo, with roughly the strength predicted by Ferrière (1996, 1998) from semi-analytic models of superbubble expansion.

All the supernova-driven models find a flow of hot gas from the plane into the halo. Avillez & Mac Low (2000) point out that much of this gas does not flow through chimneys directly into the halo, but rather rises buoyantly into the Lockman and Reynolds layers, not necessarily reaching the hot halo. A true fountain flow only appears above 1 kpc (Avillez 2000). One consequence of this is that even fairly active normal galaxies, such as the Milky Way or NGC 891, do not drive strong winds into the surrounding intergalactic medium. Kinetic energy feedback from galaxy formation will only come from starburst galaxies (Mac Low & Ferrara 1999, Silich & Tenorio-Tagle 1998, Suchkov et al. 1994, 1996).

#### Acknowledgments.

I thank the conference organizers for partial support of my attendence. This work was partially supported by the US NSF under CAREER grant AST 99-85392, and made use of the NASA Astrophysical Data System Abstract Service.

## References

Avillez, M., & Mac Low, M.-M., 2000, ApJ, in preparation

Avillez, M. 2000, MNRAS, in press (astro-ph/0001315)

Ballesteros-Paredes, J., Avillez, M. & Mac Low, M.-M. 2000, BAAS, 196, 26.08

Ballesteros-Paredes, J., Hartmann, L. & Vázquez-Semadeni, E. 1999, ApJ, 527, 285

Blitz, L. & Shu, F. H. 1980, ApJ, 238, 148

Cox, D. P. & Smith, B. W. 1974, ApJ, 189, L105

Crutcher, R. M. 1999, ApJ, 520, 706

Ferrière, K. 1996, A&A, 310, 438

Ferrière, K. 1998, A&A, 335, 488

Gerritsen, J. P. E., & Icke, V. 1997, A&A, 325, 972

Heitsch, F., Mac Low, M.-M., Klessen, R. S. 2000, ApJ, submitted

Hennebelle, P., & Pérault, M. 1999, A&A, 351, 309

Hennebelle, P., & Pérault, M. 2000, A&A, in press

Jenkins, E. B. 1978, ApJ, 219, 845

Jenkins, E. B., & Meloy, D. A. 1974, ApJ, 193, L121

Klessen, R. S., & Burkert, A. 2000a, ApJS, in press (astro-ph/9904090)

Klessen, R. S., & Burkert, A. 2000b, ApJ, in press (astro-ph/0006016)

Klessen, R. S., Burkert, A. & Bate, M. R. 1998, ApJ, 501, L205

Klessen, R. S., Heitsch, F., & Mac Low, M.-M. 2000, ApJ, 535, 887

Kim, S., Dopita, M. A., Staveley-Smith, L. & Bessell, M. S. 1999, AJ, 118, 2797

Korpi, M. 1999, PhD Thesis, University of Oulu (Finland)

Korpi, M. J., Brandenburg, A., Shukurov, A., Tuominen, I., & Nordlund, Å. 1999, ApJ, 514, L99

Mac Low, M.-M., & Ferrara, A. 1999, ApJ, 513, 142

Mac Low, M.-M., Klessen, R. S., Burkert, A., & Smith, M. D. 1998, Phys.Rev.Lett, 80, 2754

Mac Low, M.-M., & Ossenkopf, V. 2000, A&A, 353, 339

McKee, C. F. & Ostriker, J. P. 1977, ApJ, 218, 148

Passot, T., Vázquez-Semadeni, E. & Pouquet, A. 1995, ApJ, 455, 536

Rosen, A., & Bregman, J. N. 1995, ApJ, 440, 634

Scoville, N. Z., & Hersh, K. 1979, ApJ, 229, 578

Silich, S. A., & Tenorio-Tagle, G. 1998, MNRAS, 299, 249

Solomon, P. M., & Sanders, D. B. 1979, in Giant Molecular Clouds in the Galaxy, ed. P. M. Solomon & M. G. Edmunds (Oxford: Pergamon Press), 41

Stone, J. M., Ostriker, E. C. & Gammie, C. F. 1998, ApJ, 508, L99

Stutzki, J., Bensch, F., Heithausen, A., Ossenkopf, V., & Zielinsky, M. 1998, A&A, 336, 697

Suchkov, A. A., Balsara, D. S., Heckman, T. M., & Leitherner, C. 1994, ApJ, 430, 511

Suchkov, A. A., Berman, V. G., Heckman, T. M., & Balsara, D. S. 1996, ApJ, 463, 528

Nakano, T. 1998, ApJ, 494, 587

Wada, K., & Norman, C. A. 1999, ApJ, 516, L13