text:
# Large phenotype jumps in biomolecular evolution

## I Introduction

The biological function (or phenotype) of a biopolymer, such as a ribonucleic acid (RNA) or a protein, is mostly determined by the three-dimensional structure resulting from the folding of linear sequence of nucleotides (RNA) or aminoacids (proteins) that specifies a genotype. Generally, a natural biopolymer sequence (or genotype) codes for a specific two-dimensional or three-dimensional structure that defines the biopolymer activity. But one sequence can simultaneously fold in several metastable structures that can lead to different phenotypes. Thus, random mutations of a sequence induce random changes of the metastable structure populations, which generates a random walk of the biopolymer function. Understanding this phenotype random walk is a basic goal for ”quantitative” biomolecular evolution.

The statistical properties of RNA secondary structures considered as a model for genotypes have been investigated in depth in the recent years Fontana (2002). The neutral network concept Kimura (1968, 1993), i.e., the notion of a set of sequences, connected through point mutations, having roughly the same phenotype, has been shown to apply to RNA secondary structures. Thus, by drifting rapidly along the neutral network of its phenotype, a sequence may come close to another sequence with a qualitatively different phenotype, which facilitates the acquisition of new phenotypes through random evolution. Moreover, in the close vicinity of any sequence with a given structure, there exist sequences with nearly all other possible structures Schuster et al. (1994), as originally proposed in immunology Perelson and Oster (1979). Thus, even if the sequence space is much too vast to be explored through random mutations in a reasonable time (an RNA with 100 bases only has $`10^{60}`$ possible sequences), the phenotype space itself may be explored in a few mutations only, which is what matters biologically. These ideas have been brought into operation in a recent experiment Schultes and Bartel (2000) showing that a particular RNA sequence, catalyzing a given reaction, can be transformed into a sequence having a qualitatively different activity, using a small number of mutations and without ever going through inactive steps.

This paper investigates the phenotype space exploration at an elementary level by studying the statistical distribution of a population of biopolymers in a specific three-dimensional shape, that originated from a given genotypic sequence by a single mutational event. It complements studies of the evolution from one structure to another structure Fontana and Schuster (1998), that consider only the most stable structure for each sequence and neglect the thermodynamical coexistence of different structures for the same sequence. It also provides more grounds to the recent work that suggests that RNA molecules with novel phenotypes evolved from plastic populations, i.e., populations folding in several structures, of known RNA molecules Ancel and Fontana (2000). It is experimentally evident, for instance in Schultes and Bartel (2000), that some mutations change the biopolymer chemical activity by a few percents while other mutations change it by orders of magnitude. This is not unexpected since, depending on their positions in the sequence, some residues have a dramatic influence on the 3D conformation while others hardly matter. Thus, the function random walk statistically resembles a Lévy flight Shlesinger et al. (1995); Kutner et al. (1999); Bardou et al. (2002) presenting jumps at very different scales. The respective parts of gradual changes and of sudden jumps in biological evolution is a highly debated issue. While the gradualist point of view has historically dominated, evidences for the presence of jumps have accumulated at various hierarchical levels from paleontology Elredge and Gould (1972), to trophic systems, chemical reaction networks and neutral networks and molecular structure Fontana and Schuster (1998). The jump issue will be treated here by studying the statistical distribution describing the phenotype effects of random mutations of a biopolymer genotype.

To address the question of the statistical effects of random mutations of functionally active biopolymers, we propose a model inspired from disordered systems physics that naturally predicts the possibility of broad distributions of activities of randomly mutated biopolymers. With two energy parameters describing the polymer energy landscape, this models is shown to exhibit a variety of behaviors and to fit experimental data. Natural biopolymers are in a critical regime, related to the activity distribution broadness, in which a single mutation may have a large, but not too large, effect.

## II Physical model of shape population distribution

The most favorable conformational state of a biopolymer sequence with a given biological activity is generally considered to be the most stable one within the sequence energy landscape. The ruggedness of the energy landscape might vary depending on the number of other metastable, conformational states accessible by the sequence. The typical energy spacing between these states can be small enough so that several states of low energy can be populated. For simplicity, we will consider a sequence that is able to fold into its two lowest energy conformational states, an active state A of specific biological function, and an inactive state I of unknown function 111In our model, the inactive state may be replaced by an ensemble of inactive states with a given energy Onuchic et al. (1997)., but whose energy is the closest to A’s (higher or lower) (see figure 1). The differences between the free energies of the unfolded and folded states for A and I are denoted $`\mathrm{\Delta }G_\mathrm{A}`$ and $`\mathrm{\Delta }G_\mathrm{I}`$, respectively.

A mutation, i.e., a random change in the biopolymer sequence, modifies the biopolymer energy landscape so that $`\mathrm{\Delta }G_\mathrm{A}`$ and $`\mathrm{\Delta }G_\mathrm{I}`$ are transformed into $`(\mathrm{\Delta }G_\mathrm{A})_\mathrm{M}`$ and $`(\mathrm{\Delta }G_\mathrm{I})_\mathrm{M}`$. Note that the conformer state A of the mutant, its three dimensional shape, is the same as before whereas the conformer state I does not have to be the same as before. To take into account the randomness of the mutational process, the mutant free energy difference $`\delta G_\mathrm{M}(\mathrm{\Delta }G_\mathrm{I})_\mathrm{M}(\mathrm{\Delta }G_\mathrm{A})_\mathrm{M}`$ is taken either with a Gaussian distribution:

$$P_\mathrm{G}(\delta G_\mathrm{M})\frac{1}{\sqrt{2\pi }\delta G_0}e^{(\delta G_\mathrm{M}\overline{\delta G})^2/2\delta G_0^2},$$
(1)
or with a two-sided exponential (Laplace) distribution:

$$P_\mathrm{e}(\delta G_\mathrm{M})\frac{1}{2\delta G_0}e^{|\delta G_\mathrm{M}\overline{\delta G}|/\delta G_0},$$
(2)
where $`\overline{\delta G}`$ is the mean of $`\delta G_\mathrm{M}`$ and where $`\delta G_0`$ characterizes the width of the distribution. These two energy distributions are commonly used for disordered systems Doliwa and Heuer (2003) and enable us to cover a range of situations from narrow (Gaussian) to relatively broad (exponential) distributions. Assuming thermodynamic rather than kinetic control, the populations $`\pi _\mathrm{A}`$ and $`\pi _\mathrm{I}=1\pi _\mathrm{A}`$ of conformers A and I, respectively, are given by Boltzmann statistics:

$$\pi _\mathrm{A}=\frac{1}{1+e^{\delta G_\mathrm{M}/RT}},$$
(3)
where $`R`$ is the gas constant and $`T`$ is the temperature.

From the distributions of free energy differences and eq. (3), one infers the probability distributions $`P_{\mathrm{e}\mathrm{or}\mathrm{G}}(\pi _\mathrm{A})`$ of the population of conformer state A after a mutation using $`P_{\mathrm{e}\mathrm{or}\mathrm{G}}(\pi _\mathrm{A})=P_{\mathrm{e}\mathrm{or}\mathrm{G}}(\delta G_\mathrm{M})\times |\mathrm{d}\delta G_\mathrm{M}/\mathrm{d}\pi _\mathrm{A}|`$. For the Gaussian model, one obtains:

$$P_\mathrm{G}(\pi _\mathrm{A})=\frac{\mathrm{exp}\left[\frac{\left(\mathrm{ln}\pi _\mathrm{A}\mathrm{ln}(1\pi _\mathrm{A})\overline{g}\right)^2}{2g_0^2}\right]}{\sqrt{2\pi }g_0\pi _\mathrm{A}(1\pi _\mathrm{A})},$$
(4)
where $`\overline{g}\overline{\delta G}/(RT)`$, $`g_0\delta G_0/(RT)`$. The ratio $`g_0`$ of the scale of energy fluctuations and of the thermal energy appears frequently in the study of the anomalous kinetics of disordered systems. For the exponential model, one obtains:

$$P_\mathrm{e}(\pi _\mathrm{A})=\frac{e^{\overline{g}/g_0}}{2g_0\pi _\mathrm{A}^{11/g_0}(1\pi _\mathrm{A})^{1+1/g_0}}\mathrm{for}\pi _\mathrm{A}\pi _\mathrm{m},$$
(5a)
$$P_\mathrm{e}(\pi _\mathrm{A})=\frac{e^{+\overline{g}/g_0}}{2g_0\pi _\mathrm{A}^{1+1/g_0}(1\pi _\mathrm{A})^{11/g_0}}\mathrm{for}\pi _\mathrm{A}\pi _\mathrm{m},$$
(5b)
with the same definitions for $`\overline{g}`$ and $`g_0`$, and $`\pi _\mathrm{m}\left(1+e^{\overline{g}}\right)^1`$ (median population of A). Note that changing $`\overline{g}`$ into $`\overline{g}`$ is equivalent to performing a symmetry on $`P_{\mathrm{e}\mathrm{or}\mathrm{G}}(\pi _\mathrm{A})`$ by replacing $`\pi _\mathrm{A}`$ by $`1\pi _\mathrm{A}`$.

## III Types of distributions

To analyze the different types of population distributions, we focus for definiteness on the Gaussian model. A qualitatively similar behavior is obtained for the exponential model. Figure 2 represents examples of $`P_\mathrm{G}(\pi _\mathrm{A})`$ for the Gaussian model with $`\overline{g}=1`$ and various $`g_0`$’s. The negative value of $`\overline{g}`$ implies that A is on average less stable than I, and hence that $`\pi _\mathrm{A}`$ is predominantly less than 50%. For small $`g_0`$, the distribution $`P_\mathrm{G}(\pi _\mathrm{A})`$ is narrow since the width $`\delta G_0`$ of the free energy distribution is small compared to $`RT`$ so that there are only small fluctuations of population around the most probable value. When the energy broadness $`g_0`$ increases, the single narrow peak first broadens till, when $`g_01.976`$ it splits into two peaks, close respectively to $`\pi _\mathrm{A}=0`$ and to $`\pi _\mathrm{A}=1`$. The broad character of $`P_\mathrm{G}(\pi _\mathrm{A})`$ can be intuitively understood as a consequence of the non linear dependence of $`\pi _\mathrm{A}`$ on $`\delta G_\mathrm{M}`$. Thus, when the fluctuations of $`\delta G_\mathrm{M}`$ are larger than $`RT`$, i.e., when $`g_01`$, the quasi exponential dependence of $`\pi _\mathrm{A}`$ on $`\delta G_\mathrm{M}`$ (eq. (3)) non linearly magnifies $`\delta G_\mathrm{M}`$ fluctuations to yield a broad $`\pi _\mathrm{A}`$ distribution, even if $`\delta G_\mathrm{M}`$ fluctuations are relatively small compared to the mean $`\overline{\delta G}`$. A similar mechanism is at work for tunneling in disordered systems Costa et al. (2002); Romeo et al. (2003).

A global view of the possible shapes of $`P_\mathrm{G}(\pi _\mathrm{A})`$ is given in figure 3. For any given $`\overline{g}`$, when increasing $`g_0`$ starting from 0, the single narrow peak of $`P_\mathrm{G}(\pi _\mathrm{A})`$ first broadens then it splits into two peaks when $`\overline{g}=\overline{g}_{\mathrm{sign}(\overline{g})}(g_0)`$ with

$$\overline{g}_\pm (g_0)\pm \left[g_0\sqrt{g_0^22}+\mathrm{ln}\left(\frac{g_0\sqrt{g_0^22}}{g_0+\sqrt{g_0^22}}\right)\right].$$
(6)
(This expression results from a lengthy but straightforward study of $`P_\mathrm{G}(\pi _\mathrm{A})`$.) When $`g_0`$ increases further, these two peaks get closer to $`\pi _\mathrm{A}=0`$ and to $`\pi _\mathrm{A}=1`$ while acquiring significant tails (see section VII). For any given $`g_0`$, increasing $`\overline{g}`$ roughly amounts to moving the populations $`\pi _\mathrm{A}`$ towards larger values as expected since larger $`\overline{g}`$’s correspond to stabler states A. However, distinct behaviours arise depending on $`g_0`$. If $`g_0<\sqrt{2}`$, whatever the value of $`\overline{g}`$, the distribution $`P_\mathrm{G}(\pi _\mathrm{A})`$ is always sufficiently narrow to present a single peak. If $`g_0>\sqrt{2}`$, the distribution $`P_\mathrm{G}(\pi _\mathrm{A})`$ is sufficiently broad to have two peaks when, furthermore, the distribution is not too asymmetric, which occurs for $`\overline{g}[\overline{g}_{}(g_0),\overline{g}_+(g_0)]`$. In short, depending on $`\overline{g}=\delta \overline{G}/RT`$, which characterizes mainly the peak(s) position, and on $`g_0=G_0/RT`$, which characterizes mainly the distribution broadness, the distributions $`P_\mathrm{G}(\pi _\mathrm{A})`$ are either unimodal or bimodal, either broad or narrow. This variety of behaviors is reminiscent of beta distributions.

## IV From shape populations to catalytic activities

Up to now, we have discussed the distribution $`P(\pi _\mathrm{A})`$ of the population of a shape A that is functionally active. However, as far as it concerns biopolymers with enzymatic functions, what is usually measured is a chemical activity $`a`$, i.e., the product of a reaction rate $`k`$ for the conformer A by the population $`\pi _\mathrm{A}`$ of this conformer. The reaction rates are given by the Arrhenius law $`k=k_0e^{E_\mathrm{a}/RT}`$ where $`k_0`$ is a constant and $`E_\mathrm{a}`$ is the activation energy. Thus, the chemical activity writes, using eq. (3):

$$a=k_0e^{E_\mathrm{a}/RT}\frac{1}{1+e^{\delta G_\mathrm{M}/RT}}.$$
(7)
Random mutations may induce random modifications of $`E_\mathrm{a}`$, $`\delta G_\mathrm{M}`$ or both. Fluctuations of $`\delta G_\mathrm{M}`$ have been treated above. One can introduce fluctuations of $`E_\mathrm{a}`$ in the same way. We do not do it here in details but present only the general trends.

The effects of adding an activation energy distribution in addition to the free energy difference distribution are twofold. For small activities, the distribution $`P(a)`$ of chemical activities is similar to the small $`P(\pi _\mathrm{A})`$ peak at small $`\pi _\mathrm{A}`$. Indeed, the reaction rate $`k`$ depends exponentially on $`E_\mathrm{a}`$, just as the population $`\pi _\mathrm{A}`$ depends exponentially on $`\delta G_\mathrm{M}`$ when $`\pi _\mathrm{A}1`$. Moreover, the product of two broadly distributed random variables is also broadly distributed 222With Gaussian distributions of $`\delta G_\mathrm{M}`$ and $`E_\mathrm{a}`$, one can be more specific. Both $`\pi _\mathrm{A}`$ and $`k`$ are then lognormally distributed at small values. Thus, the product $`a=k\pi _\mathrm{A}`$ is also lognormally distributed Romeo et al. (2003). with a shape similar to the one of $`P(\pi _\mathrm{A})`$. For large activities, on the other hand, $`\pi _\mathrm{A}`$ and $`k`$ behave differently because $`\pi _\mathrm{A}`$ is bounded by 1 while $`k`$ is unbounded. Thus, if the $`k`$ distribution is broad enough, the distribution of $`a`$ at large $`a`$ may exhibit a broadened structure compared to the $`\pi _\mathrm{A}1`$ peak of $`P(\pi _\mathrm{A})`$.

In summary, the distribution of chemical activities $`P(a)`$ is similar to the distribution of shape populations $`P(\pi _\mathrm{A})`$ when $`P(\pi _\mathrm{A})`$ presents a large $`\pi _\mathrm{A}0`$ peak (conditions for this to occur are explicited in section VI). Thus, by observing the shape of the $`a0`$ peak in the activity distribution $`P(a)`$, one does not easily distinguish between activation energy dispersion, which affects $`k`$, and free energy difference dispersion, which affects $`\pi _\mathrm{A}`$. On the other hand, at large $`a`$, $`P(a)`$ is differently influenced by activation energy dispersion and by free energy difference dispersion. The available experimental data (see section V) enables us to analyze precisely $`P(a)`$ at small activities but not at large activities. Thus, for practical purposes, it is not meaningful in this paper to consider a distribution of activation energies on top of a distribution of free energy differences. In the sequel, we will thus do as if only the distribution of free energies was involved, stressing that similar effects can be obtained from a distribution of activation energies.

## V Analysis of experimental data

Comparison of the theoretical distributions of eq. (5) and eq. (4) with experimental data enables us to test the relevance of the proposed model. We have analyzed the measurements of the catalytic activities of a set of 157 mutants derived from a self-splicing group I ribozyme, a catalytic RNA molecule Couture et al. (1990) (out of the 345 mutants generated in Couture et al. (1990), we only considered the 157 ones with single point mutations). The original ”wild-type” molecule is formed of a conserved catalytic core that catalyzes the cleavage of another part of the molecule considered as the substrate. The set of mutants is derived from the original ribozyme by performing systematically all single point mutations of the catalytic core, i.e., of the part of the molecule that influences most the catalytic activity. Nucleotides out of the core, that in general influence less the catalytic activity, are left unmutated. Thus, in our framework, this set of mutants can be seen as biased towards deleterious mutations. Indeed, mutations of the quasi optimized core are likely to lead to much less active mutants, while mutations of remote parts are likely to leave the activity essentially unchanged. If all parts of the molecule had been mutated, more neutral or quasi neutral mutations would have been obtained. Another point of view, which we adopt here is to consider the catalytic core as a molecule in itself, on which all possible single point mutations have been performed.

The 157 measured activities are used to calculate a population distribution with inhomogenous binning (cf. broad distribution). Two bins required special treament: the smallest bin, centered in 0.5%, contains 40 mutants with non measurably small activities ($`<1\%`$ of the original activity); the largest bin, centered in 95%, contains the 6 mutants with activities larger than 90% of the original ’wild’ RNA activity (the largest measured mutant activity is 140%). These two points, whose abscissae are arbitrary within an interval, are not essential for the obtained results. At last, as very few mutants have activities larger than the wild-type ribozyme, the proportionality constant between activity and population is set by matching a population $`\pi _\mathrm{A}=1`$ to the activity of the wild-type ribozyme.

The obtained distribution (see figure 4) has a large peak in $`\pi _\mathrm{A}0`$, indicating that most mutations are deleterious, with a long tail at larger activities and a possible smaller peak in $`\pi _\mathrm{A}1`$. This non trivial shape is well fitted by the Gaussian model of eq. (4) with $`\overline{g}=3.6`$ and $`g_0=2.9`$ (the uncertainty on these parameters is about 50%, see dashed lines in figure 4). One infers $`\overline{\delta G}2.1`$ kcal/mol and $`\delta G_01.7`$ kcal/mol ($`T=300`$ K). The order of magnitude of these values is compatible with thermodynamic measurements performed on similar systems Jaeger et al. (1993, 1994); Brion and Westhof (1997); Brion et al. (1999). This confirms the plausibility of the proposed approach. The inset of figure 4 shows the population distribution in the exponential model with $`\overline{g}`$ and $`g_0`$ values taken from the Gaussian fit. The agreement with the experimental data is also quite good. Thus, the proposed approach soundly does not strongly depend on the yet unknown shape details of the energy distribution. Finally, one can estimate the broad character of the activity distribution from the statistical analysis of the experimental data. Indeed, according, e.g., to the Gaussian model fit, the typical, most probable, population $`\pi _\mathrm{A}`$ is found to be $`6\times 10^6`$ while the mean population is $`0.15`$. Thus, the activity distribution spans more than four orders of magnitude.

## VI Coarse graining description: all or none features

The variation of activity of a biopolymer upon mutation is often described as an ‘all or none’ process: mutations are considered either as neutral (the mutant retains fully its activity and $`\pi _\mathrm{A}100`$ %) or as lethal (the mutant loses completely its activity and $`\pi _\mathrm{A}0`$ %). Satisfactorily, a coarse graining description of the proposed statistical models exhibits such all or none regimes for appropriate $`(\overline{g},g_0)`$ values, as well as other regimes.

To obtain a quantitative coarse graining description, we define the mutants with ’no’ activity as those with population that has less than $`12\%`$ ($`\pi _\mathrm{A}(\delta G_\mathrm{M}=2RT)`$) in the A shape. Their weight is

$$w_0=_0^{12\%}P_{\mathrm{e}\mathrm{or}\mathrm{G}}(\pi _\mathrm{A})d\pi _\mathrm{A}=_{\mathrm{}}^{2RT}P_{\mathrm{e}\mathrm{or}\mathrm{G}}(\delta G)d\delta G.$$
(8)
Similarly, the mutants with ‘full’, respectively ’intermediate’, activity are defined as those with $`\pi _\mathrm{A}88\%`$, respectively $`12\%\pi _\mathrm{A}88\%`$, and their weight is $`w_{100}=_{2RT}^{\mathrm{}}P_{\mathrm{e}\mathrm{or}\mathrm{G}}(\delta G)d\delta G`$, respectively $`w_\mathrm{i}=_{2RT}^{+2RT}P_{\mathrm{e}\mathrm{or}\mathrm{G}}(\delta G)d\delta G`$. Taking for definiteness the Gaussian model leads to

$$w_0=\mathrm{\Phi }\left(\frac{2}{g_0}\frac{\overline{g}}{g_0}\right)$$
(9)
where $`\mathrm{\Phi }(u)=_{\mathrm{}}^ue^{t^2/2}dt/\sqrt{2\pi }`$ is the distribution function of the normal distribution. Similarly, one has $`w_\mathrm{i}=\mathrm{\Phi }[(2\overline{g})/g_0]\mathrm{\Phi }[(2+\overline{g})/g_0]`$ and $`w_{100}=1\mathrm{\Phi }[(2\overline{g})/g_0]`$. Approximate expressions for $`\mathrm{\Phi }(u)`$ ($`\mathrm{\Phi }(u)e^{u^2/2}/(\sqrt{2\pi }u)`$ for $`u1`$, $`\mathrm{\Phi }(u)1/2+u/\sqrt{2\pi }`$ for $`|u|1`$ and $`\mathrm{\Phi }(u)1e^{u^2/2}/(\sqrt{2\pi }u)`$ for $`u1`$) give the regimes in which each weight $`w`$ is negligible ($`w1`$), dominant ($`1w1`$) or in between. For instance, $`w_0`$ is negligible for $`\overline{g}>g_02`$, dominant for $`\overline{g}<g_02`$ and intermediate for $`g_02<\overline{g}<g_02`$. These inequalities indicate the transition from one regime to another. To be strictly in one regime requires typically that $`\overline{g}/g_0`$ is larger or greater than 1 from the corresponding criterion, e.g., $`w_0`$ is strictly negligible when $`\overline{g}/g_0>1+(g_02)/g_0`$ The transitions from one regime to another one are in general exponentially fast (solid lines in figure 5). However, in the region ($`g_0>2`$, $`|\overline{g}|<g_02`$), the transitions from one regime to another one are smooth (dashed lines in figure 5) since, in this region, the weights vary slowly, e.g., $`w_\mathrm{i}4/(g_0\sqrt{2\pi })`$.

The resulting coarse graining classification of $`P_\mathrm{G}(\pi _\mathrm{A})`$ is represented in figure 5. The ‘all or none’ behaviour, denoted ‘0 & 100’, appears in the region $`g_06/\sqrt{\pi /2}`$ and $`|\overline{g}|\sqrt{\pi /2}g_06`$ as the result of a large dispersion of energy differences associated to a moderate average energy difference. We note that all possible types of distributions are actually present in this model: probabilities concentrated at small, intermediate or large values (0, i or 100); probabilities spread over both small and intermediate (0 & i), both small and large (0 & 100, all or none) or both intermediate and large (i & 100) values; probabilities spread over small, intermediate and large values at the same time (0 & i & 100). The coarse graining classification of figure 5 complements the number of peaks classification of figure 3 without overlapping it. Indeed, there exist parameters $`g_0`$ and $`\overline{g}`$ for which, e.g., two peaks coexist but one of these peaks has a negligible weight. Thus the presence of a peak is not automatically associated to a large weight in the region of this peak.

## VII Zooming in the $`\pi _\mathrm{A}0`$ peak: long tails

To go beyond the coarse graining description, we zoom in the $`\pi _\mathrm{A}0`$ peak. As shown in the inset of figure 3, the small activities, labelled as ‘no activity’ in a coarse graining description, actually consist of non zero activities with values scanning several orders of magnitude. This can be analyzed quantitatively, e.g., in the Gaussian model. For $`\pi _\mathrm{A}0`$, the activity distribution given by eq. (4) is quasi lognormal:

$$P_\mathrm{G}(\pi _\mathrm{A})\frac{1}{\sqrt{2\pi }g_0\pi _\mathrm{A}}\mathrm{exp}\left[\frac{(\mathrm{ln}\pi _\mathrm{A}\overline{g})^2}{2g_0^2}\right].$$
(10)
Thus, $`P_\mathrm{G}(\pi _\mathrm{A})`$ has as a power law like behavior Romeo et al. (2003); Montroll and Shlesinger (1983):

$$P_\mathrm{G}(\pi _\mathrm{A})\frac{1}{\sqrt{2\pi }g_0\pi _\mathrm{A}}\mathrm{for}e^{\overline{g}\sqrt{2}g_0}\pi _\mathrm{A}e^{\overline{g}+\sqrt{2}g_0},$$
(11)
in the vicinity of the lognormal median $`e^{\overline{g}}`$. This corresponds to an extremely long tailed distribution, since $`1/\pi _\mathrm{A}`$ is not even normalizable. It presents the peculiarity that, for $`a`$ and $`a+1`$ belonging to $`[\overline{g}\sqrt{2}g_0,\overline{g}+\sqrt{2}g_0]`$, the probability to obtain a population $`\pi _\mathrm{A}`$ of a given order of magnitude $`a`$, i.e., $`\pi _\mathrm{A}[e^a,e^{a+1}]`$, does not depend on the considered ordered of magnitude $`a`$, since

$$_{e^a}^{e^{a+1}}P_\mathrm{G}(\pi _\mathrm{A})d\pi _\mathrm{A}\mathrm{const}.$$
(12)
Thus, if a living organism has to adapt the chemical activity of one of its biopolymer constituents, it can explore several order of magnitude of activity by only few mutations within the biopolymer. The activity changes mimic a Lévy flight Bouchaud and Georges (1990) as revealed, e.g., by the experimental data in Schultes and Bartel (2000). The large activity changes will raise self-averaging issues Romeo et al. (2003) that will add up to those generated by correlations along evolutionary paths Bastolla et al. (2002)

Three broadness regimes corresponding to three evolutionary regimes can be distinguished. If $`g_0`$ is very large, the mutant activities span a very large range. This regime might be globally lethal because, in most cases, the mutant activity will be either too low or too large to be biologically useful. However, under conditions of intense stress, the large variability might allow the system to evolve radically. With $`g_0=10`$, for instance, the activity range covers typically 12 orders of magnitude from $`10^6e^{\overline{g}}`$ to $`10^6e^{\overline{g}}`$ (see eq. (11)). If $`g_0`$ is moderately large, the mutant activities span just a few orders of magnitude. This regime is broad enough to permit significant changes, but not too broad to avoid producing too many lethal changes. With $`g_0=3`$, for instance, the activity range covers typically $`34`$ orders of magnitude from $`10^{1.8}e^{\overline{g}}`$ to $`10^{1.8}e^{\overline{g}}`$. If $`g_0`$ is small, the lognormal distribution peak can be approximated by a Gaussian Romeo et al. (2003)

$$P_\mathrm{G}(\pi _\mathrm{A})\frac{1}{\sqrt{2\pi }g_0e^{\overline{g}}}\mathrm{exp}\left[\frac{(\pi _\mathrm{A}e^{\overline{g}})^2}{2(g_0e^{\overline{g}})^2}\right].$$
(13)
The distribution is now narrow and the ranges of values is typically $`[e^{\overline{g}}(12g_0),e^{\overline{g}}(1+2g_0)]`$. This type of distribution is not adapted for producing large changes, but rather for performing fine tuning optimization. With $`g_0=0.1`$, for instance, the activity range covers only $`\pm 20\%`$ around $`e^{\overline{g}}`$.

We remark that the group I ribozyme which we have analyzed corresponds to $`g_02.9`$, right in the critical regime of moderately large $`g_0`$. One can guess from experimental studies of other biopolymers or from chemical considerations that most biopolymers will fall in this range since $`\delta G_0`$ is typically on the order of a few kilocalories while $`RT`$ is $`0.6`$ kcal (Note that $`\delta G_0`$ corresponds to the free energy change between the biopolymer native 3D state and an unfolded state, in which the biopolymer has lost its three-dimensional shape but not its full secondary structure). It would be interesting to perform further statistical data analysis to see how, e.g., the available protein mutagenesis studies fit with our present model.

The energy statistics associated mutations is likely to be determined at gross scale by the basic biophysics of the molecules involved. This fixes a range for $`\delta G_0`$. It is nonetheless plausible and suggested by our discussion that there is an evolutionary preferred type of activity distribution, and hence of sequences, that may imply a fine tuning of $`g_0=\delta G_0/RT`$ within the constraints on $`\delta G_0`$ coming from biophysics (see figure 6) so that each mutation typically generates a significant, but not systematically lethal, activity change. If one considers that the activity changes must cover between, say, one and seven orders of magnitude, then the allowed $`g_0`$ range is $`16`$ (see eq. (11)).

To answer the question whether the energy statistics is solely dictated by molecular biophysics or whether it is also influenced by evolutionary requirements, one may compare the energy statistics of molecules from different thermal environments. The conservation of the $`\delta G_0`$ range across psychrophilic and thermophilic molecules would stress the domination of biophysics factors. Note that our model would then imply different stochastic evolutionary dynamics, through the width of the activity distribution, for psychrophilic and thermophilic environments. Conversely, the conservation of $`g_0`$ would reveal the importance of evolutionary requirements.

## VIII Conclusions

In this paper, we have presented a model for the distribution of biopolymer activities resulting from mutations of a given sequence. The model is characterized by the statistics of the energy differences between active conformations and inactive conformations. A similar model would be obtained by considering the statistics of activation energies. The model fits the measured activity distribution of a ribozyme with energy parameters in the physically appropriate range. It is also able to reproduce commonly observed behaviours such as all or none.

Importantly, the peak of small activities exhibits three distinct types depending on the broadness of the distribution of energy differences. Real biopolymers are in a critical regime allowing the exploration of different ranges of activities in a few mutations without being too often lethal. This critical regime seems the most favorable evolutionary regime and could be the statistical engine allowing molecular evolution. Thus the present work supports the idea that, for evolution to take place, the temperature and the physico-chemistry dictating the free energy scales of biopolymers must obey a certain ratio. At last, it suggests that, by looking at small variations of this ratio, one might be able to classify biopolymers. One expects, for instance, that biopolymer sequences that are locked in a shape with a specific function, will have smaller $`g_0`$ than rapidly evolving biopolymers sequences that could acquire new functions by undergoing major structural changes. Thus, at the origin of life or during rapidly evolving punctuations, biopolymers with larger $`g_0`$ than those characterizing highly optimized, modern RNA and protein molecules, could have contributed to the emergence of novel phenotypes, leading thus to an increase of complexity.


format: text

text:
# The first detection of Far-Infrared emission associated with an extended HI diskBased on observations with the Infrared Space Observatory (ISO), an ESA project with instruments funded by ESA member states (especially the PI countries: France, Germany, the Netherlands, and the United Kingdom) and with the participation of ISAS and NASA.. The case of NGC 891.

## 1 Introduction

Radial gas surface density profiles in spiral galaxies show quite similar behaviour in relation to the optical disk, irrespective of their morphological type (Sancisi 1995, 1999). Whereas the molecular gas is concentrated towards the inner disk, the HI surface density is generally flat at an average level of $`10\mathrm{M}_{}\mathrm{pc}^2`$ over the entire extent of the optical disk (though with some variation from galaxy to galaxy, see Broeils & van Woerden 1994). Exterior to the optical disk, which tends to have a comparatively abrupt cut-off at $`3`$ stellar exponential scale lengths (Pohlen et al. 2000), the HI surface density falls off exponentially until a level of ca. $`0.1\mathrm{M}_{}\mathrm{pc}^2`$ is reached. At this point the gas disk either ends, or becomes ionised by the intergalactic medium. The portion of the HI disk extending beyond the optical stellar disk is commonly referred to as the “extended HI disk”.

It is unknown whether these gaseous disks are remnants of primordial material left over from the epoch of galaxy formation (Larson 1990), or whether they contain material reprocessed in stellar interiors, either transferred from the stellar disk or captured from other galaxies. In the latter case the extended HI disks should be enriched by metals produced in stars (Tinsley & Larson 1978, Pei, Fall & Hauser 1999, Maller et al. 2001), so observations of these species could be used to identify their nature. Unfortunately metals in extended HI disks are difficult to detect in the gas phase, either through their emission line spectrum (because of the lack of exciting stars) or through absorption (because of the lack of sufficient background sources). However any metals present in form of dust grains offers an alternative way to trace the origin of these gaseous disks. Tentative evidence for the presence of metals in the form of grains was provided from measurements of colour variations between background galaxies and control fields (Zaritsky 1994). There is also evidence that, within the confines of the optical disk, grains have a larger scale length than the stars (Alton et al. 1998a, Davies et al. 1999, Trewhella et al. 2000, Radovich, Kahanpää & Lemke 2001, Xilouris et al. 1998, Xilouris et al. 1999, Popescu et al. 2000a, Misiriotis et al. 2001) and that dust extends right up to the edge of the optical disks (Cuillandre et al. 2001). In some Blue Compact Dwarf galaxies observed in the FIR by Tuffs et al. (2002a,b) it has also been suggested that there is dust outside the optical emitting core region (Popescu et al. 2002).

Here we present the first detection of cold dust in an extended HI disk, achieved through a dedicated deep Far-Infrared (FIR) observation of a large field encompassing the entire HI disk of the edge-on spiral galaxy NGC 891, made using the ISOPHOT instrument (Lemke et al. 1996) on board the Infrared Space Observatory (ISO)(Kessler et al. 1996). We chose NGC 891 for this observation as it has an asymmetric HI disk (Swaters et al. 1997), so any FIR counterpart should also be asymmetric and thus be more easily recognisable.

## 2 Observations and data reduction

A detailed description of the observations and data reduction is given in Popescu et al. (2003). Here we give only a brief overview of these procedures. The observations were made using the ISOPHOT-C200 2$`\times 2`$ pixel array in the C160 and C200 filters, which respectively cover passbands of $`130\mathrm{\hspace{0.17em}218}`$ and $`170\mathrm{\hspace{0.17em}239}\mu \mathrm{m}`$ and have central wavelengths of 170 and 200 $`\mu `$m. These FIR wavelengths were chosen since they provide maximum sensitivity to any cold dust present in the extended HI disk. Since there are no local heating sources in this region and the only available photons to heat the grains would be those coming from the optical disk, the FIR emission from any embedded dust was anticipated to be faint and have a spectral peak near our chosen filters. In order to cover the entire HI disk of NGC 891 (extending up to $`10`$ arcmin from the nucleus; Swaters et al. 1997), as well as the surrounding background, a field of radius $`\pm 13.5`$ arcmin ($`\pm 40`$ kpc) was mapped along the major axis of the galaxy.

The FIR maps were obtained using the “P32” mapping mode which provided near Nyquist sampling over the three overlapping fields: north, south and central. Preliminary results obtained on the central field were presented by Popescu & Tuffs (2002a), Tuffs & Popescu (2003) and by Dupac et al. (2003). The data were processed using the latest P32 reduction package (Tuffs & Gabriel 2003), which corrects for the transient response of the detector pixels. This allowed high dynamic range maps to be constructed to levels of 1 percent of the peak disk brightness. A time dependent flat field correction was made for each map, by fitting a cubic function to the response of the detector pixels to the background. Calibration was made using V8.1 of the ISOPHOT Interactive Analysis Package PIA (Gabriel et al. 1997).

Radial profiles were derived from the background subtracted FIR maps by integrating the emission parallel to the minor axis of the galaxy in bins of width 18 arcsec along the major axis. Independent data contribute to each map pixel (and to each point on the profile)(Tuffs & Gabriel 2003)

## 3 Results

NGC 891 was found to be a normal galaxy in respect to its integrated FIR properties, such that the fraction of stellar light reradiated by dust in the FIR is $`30\%`$ (Popescu & Tuffs 2002a), which is close to the mean value for normal galaxies (Popescu & Tuffs 2002b). A detailed analysis of the FIR surface brightness distribution and profiles within the optical disk ($`\pm 360^{\prime \prime }`$) has been presented in Popescu et al. (2003). There we find that the measured profiles and surface brightness distributions are in excellent agreement with the prediction for their counterparts obtained using the model for the optical/FIR/submm spectral energy distribution of Popescu et al. (2000a). The derived intrinsic distributions of dust and stars in NGC 891 were constrained from the optical/NIR images of NGC 891, as well as using data at 60, 100, 450 & 850 $`\mu `$m (Alton et al. 1998b) and data at 1.3 mm (Guelin et al. 1993). Here we compare the FIR profiles with the HI profiles and concentrate on the emission beyond the optical disk.

In Fig. 1 we show the resulting radial profile at 200 $`\mu `$m, overlaid with the corresponding HI profile. The latter was obtained from the HI maps of Swaters et al. (1997) after convolution to the PSF of the ISO measurements. Within $`\pm 200^{\prime \prime }`$ from the centre, where the HI radial profile is fairly flat, the 200 $`\mu `$m profile rises continuously towards the nucleus. This can be attributed in part to an increasing surface density of grains associated with molecular gas, which is known to predominate in the inner disk (García-Burillo & Guélin 1995), as well as to the stronger radiation fields from the inner parts of the galaxy. Between $`\pm 200^{\prime \prime }`$ and $`\pm 360^{\prime \prime }`$ (the edge of the optical disk) both the FIR and the HI profiles fall steeply. The $`200\mu `$m profile can be traced out to a radius of $`522^{\prime \prime }`$ (24 kpc) in the South - $`160^{\prime \prime }`$ (7.4 kpc) beyond the edge of the optical disk, and out to $`432^{\prime \prime }`$ (19.9 kpc) in the North - $`70^{\prime \prime }`$ beyond the optical disk. By comparison the HI profile extends out to $`600^{\prime \prime }`$ in the South and out to $`440^{\prime \prime }`$ in the North (Swaters et al. 1997). Thus, the extent and the asymmetry in the 200 $`\mu `$m emission follows that of the HI emission (the same is true for the 170 $`\mu `$m radial profile), indicating a dust emission counterpart to the extended HI disk.

In order to check that the FIR emission detected in the extended HI disk is not attributable to beam smearing of emission from within the outer confines of the optical disk, we truncated the HI distribution at the edge of the optical disk (such that no emission would exist beyond this edge), convolved it with the ISOPHOT beam, and overplotted the result in Fig. 1. It is clear that the FIR emission detected (in the South) beyond the optical disk falls above this modified HI profile, indicating that the detection is real and not of an instrumental nature. We also note that the bright emission from the central region of the galaxy would not contribute to the emission seen towards the extended HI disk, as the ISOPHOT beam is known to fall steeply, and has no extended wings.

The relative rate at which the 200 $`\mu `$m and HI profiles decline with radius is best seen in Fig. 2 (top panel), where the profile of the ratio of these quantities is plotted. On the southern side this ratio decreases by a factor of 2.3 between $`360^{\prime \prime }`$ (the edge of the optical disk) and $`522^{\prime \prime }`$ (the maximum extent of the detected FIR emission). This could either be due to a decrease in the dust-to-gas ratio, or to a decrease in the grain heating with radius, or to a combination of the two. The dust-to-gas ratio can be calculated from the ratio plotted in Fig. 2 (top panel), if the dust temperature $`T_\mathrm{D}`$ is known. In Fig. 2 (bottom panel) it is shown that the 200/170 radial colour profile (measured as far as a radius of $`400^{\prime \prime }`$) has a smooth progression towards colder emission with increasing radial distance, which can be taken as evidence that the dust at large galactic radii should be cold. From the 200/170 colour ratio of 1.1, as measured at a radius of $`400^{\prime \prime }`$, and assuming a grain emissivity proportional to $`\lambda ^2`$, we derive $`T_\mathrm{D}=14`$ K. This provides us with an upper limit for $`T_\mathrm{D}`$ at radii larger than $`400^{\prime \prime }`$, because the radiation fields should decrease continuously with increasing distance from the edge of the optical disk. For the distance of NGC 891 of 9.5 Mpc (van der Kruit & Searle 1981) and assuming a mixture of silicate and graphite grains (Laor & Draine 1993, Draine & Lee 1984), we derive a dust-to-gas mass ratio of 0.0083 at the projected radius of $`400^{\prime \prime }`$. This value is close to the dust-to-gas ratio of $`1\%`$ of both the local Galactic interstellar medium and of NGC 891 as a whole. Thus our measurements are consistent with there being little or no variation in the dust-to-gas ratio between the optical disk and the extended HI disk. So the radial fall-off in the ratio between the 200 $`\mu `$m and HI profiles in the extended HI disk is most probably due to the diminishing heating rate. This inference falls broadly in line with the expected decline in photon flux by a factor 2.1, if the photons heating the grains originate from near the edge of the optical disk.

## 4 Discussion and Conclusions

The existence of large amounts of grains in the extended HI disk of NGC 891 raises the challenging question about their origin and the implications for the origin of the extended HI disk itself. The large value of the dust-to-gas ratio obtained for the extended HI disk clearly indicates that this gaseous disk is not primordial, left over from the epoch of galaxy formation. The detected grains must have either been transported from the optical disk, or they must have been produced outside the galaxy.

If the grains were transported from the optical disk continuously over the lifetime of the galaxy, only a very small fraction of grains produced in the optical disk need to be transferred to explain the derived dust mass in the extended HI disk, since there are no obvious grain destructions mechanisms operating there. Taking the lifetime of NGC 891 to be $`\tau _{\mathrm{gal}}=10^{10}`$ yr and the survival timescale of grains in the optical disk $`\tau _{\mathrm{surv}}=10^8`$ yr, the total grain production mass is $`M_D^{\mathrm{tot}}=(\tau _{\mathrm{gal}}/\tau _{\mathrm{surv}})\times M_\mathrm{D}^{\mathrm{opt}}`$, where $`M_\mathrm{D}^{\mathrm{opt}}`$ is the observed mass of grains in the optical disk at the current epoch (Popescu et al 2000a). We obtain $`M_D^{\mathrm{tot}}=1.3\times 10^{10}M_{}`$. If we compare this mass with the total mass of dust present in the extended HI disk $`M_D^{\mathrm{ext}}=2.3\times 10^6\mathrm{M}_{}`$ we can derive an efficiency of transfer of dust grains from the optical disk $`\eta =1.8\times 10^4`$. Thus, it would only require a tiny amount of grains to be transported to the extended HI disk to account for the observations. This raises the possibility that the grains were transported via the prominent halo of NGC 891, as originally traced in H$`\alpha `$ by Dettmar (1990) and Rand et al. (1990). For the specific case of transporting grains into the halo several mechanisms have been proposed by Ferrara (1991), Davies et al. (1998) and Popescu et al. (2000b), though no theory exists for the transport of grains through the halo to higher galactocentric radii. However it would be a remarkable coincidence that the transfer efficiency inferred from our observations should take exactly the value for which the present dust-to-gas ratio in the extended HI disk matches the dust-to-gas ratio in the optical disk. Furthermore any transport of grains via the halo should produce a symmetrical distribution of dust, contrary with what is observed.

An alternative mechanism for transporting grains and gas from the inner disk would be diffusion triggered by macro turbulence. To explain our observation, this mechanism would also have to be effective beyond the optical disk, in regions unperturbed by mechanical energy input from supernovae and stellar winds. A further requirement for this mechanism to be effective would be that the timescale for the mixing mechanism should be shorter than the timescale for grain destruction in the optical disk. As in the case of the transport of grains via the halo, an argument against this mechanism is however the observed asymmetry of the FIR profile.

Another possibility is that both the gas and the dust in the extended HI disk were part of the interstellar medium of another galaxy which was (long ago) tidally stripped and captured by NGC 891111We note that NGC 891 is thought to be a non-interacting system at the current epoch.. This would also explain the asymmetry in both the HI and in the dust. A present day example of such an interaction-accretion event is the advanced interaction of a dwarf galaxy with M 101 (van der Hulst & Sancisi 1988).

To conclude, while the exact mechanism through which the extended HI disk is formed remains unclear, our detection of FIR emission rules out a primordial origin of the extended HI disk in NGC 891. For the moment our result was obtained for one galaxy and cannot therefore be generalised to all spiral galaxies. However, future observations will be able to prove if in general extended HI disks of spiral galaxies contain large amounts of dust or if this is a characteristic peculiar to NGC 891. If the former is the case, then this implies that at the epoch when the first galaxies formed, there must have been a rather efficient process which removed primordial debris from around the forming galaxies, for example in a strong galactic wind, or simply as a result of a rather efficient conversion of gas into stars. Another implication of our detection of an asymmetric dust counterpart to the extended HI disk in NGC 891 is that the asymmetry of the latter is intrinsic rather than being due to the disk becoming ionised at a shorter radius.

###### Acknowledgements.

We would like to thank the anonymous referee for his useful comments and suggestions. We would also like to thank Andreas Burkert and Jay Gallagher for helpful discussions and Rob Swaters for providing us with the HI map.


text:
# A Resolved Circumstellar Disk around the Herbig Ae Star HD 100546 in the Thermal Infrared

## 1 Introduction

Circumstellar disks provide insight into the formation of planetary systems. These disks are observed most readily around luminous pre-main-sequence (PMS) stars. Herbig Ae (HAE) stars, the evolutionary precursors to intermediate mass main-sequence stars such as Vega, have been identified to have infrared (IR) excess emission. The source of the emission has been hypothesized by Hillenbrand et al. (1992) and Lada & Adams (1992) to originate from a geometrically thin, optically thick circumstellar disk, with an optically thin inner region and a high accretion rate in order to explain the observed spectral energy distribution (SED) of such stars. An alternative interpretation suggests that the emission may be a result of a ”dusty nebula” (or envelope), rather than a disk (Hartmann, Kenyon, & Calvet, 1993). Recent modelling has shown a likely possibility to be a disk which flares vertically with increasing radius from the star (Chiang & Goldreich, 1997; Kenyon & Hartmann, 1987). Other studies have found that these disks may have a more complex structure, incorporating an inner hole and heating of the inner wall of the disk to account for an excess in near-IR emission (Dullemond, Dominik, & Natta, 2001). Another model incorporates an extended spherical envelope surrounding a thin disk (Miroshnichenko et al., 1999). We refer the reader to Natta, Grinin, & Mannings (2000) for an extensive review of recent results. In general, observations of the circumstellar environments of PMS stars are an important step in determining which models are most representative of their true environment, as well as understanding the evolution of protoplanetary disks into planetary systems.

The nearby ($``$100 pc) HAE star HD 100546 has been the focus of several studies. Malfait et al. (1998) characterized the spectrum of the star in the IR, and identified several spectral features indicative of silicate and polycyclic aromatic hydrocarbon (PAH) species in the circumstellar environment. They also found features in the spectrum of HD 100546 to be very similar to those in comet Hale-Bopp, indicating the presence of cometary material in the system, and hypothesize that the system could harbor giant protoplanets to explain the presence of crystalline silicates in the cometary material. A recent study by Bouwman et al. (2003) found the spectrum of HD 100546 to be dramatically different from other HAE stars, and propose a model with a circumstellar disk with an inner gap of 10 AU and a giant protoplanet. Three studies, Grady et al. (2001); Augereau et al. (2001) and Pantin, Waelkens, & Lagage (2000), used coronagraphic observations at near-IR wavelengths to image the dust disk in scattered light and characterize its spatial structure. These studies detect evidence of an inclined dust disk, and are in good agreement as to its inclination ($`40\mathrm{°}`$ from face-on), and position angle of its semimajor axis ($`130\mathrm{°}`$ to $`160\mathrm{°}`$ E of N). Extended emission has also been detected at 3.4 mm (Wilner et al., 2003), and far-ultraviolet observations of warm molecular hydrogen are also consistent with the presence of an inclined disk (Lecavelier des Etangs et al., 2003). The presence of circumstellar emission from HD 100546, as well as its relative proximity, make it an ideal target for nulling interferometry.

Nulling interferometry is a technique used to study circumstellar environments by suppressing the starlight which normally overwhelms any signal from the faint circumstellar material. The technique is implemented by overlapping the pupils of two telescopes (or two subapertures from a single telescope) with an appropriate path length difference to destructively interfere the light. The result is a sinusoidal transmission pattern where the unresolved central point source is suppressed and the surrounding resolved structure can be detected.

In this Letter, we present the results of nulling interferometric and direct imaging observations of the HAE star HD 100546 in the mid-IR. We discuss the detection and structure of resolved emission surrounding the star at several wavelengths. For further background we refer the reader to Hinz, Hoffmann, & Hora (2001).

## 2 Observations & Data Reduction

Observations were made in August 2001 and May 2002 at the Magellan I (Baade) 6.5 m telescope at Las Campanas Observatory near La Serena, Chile. The BracewelL Infrared Nulling Cryostat (BLINC) was the nulling interferometer used to provide suppression of starlight, creating an interferometer with two 2.5 m subapertures and a baseline of 4 m. The Mid-InfraRed Array Camera (Hoffmann et al., 1998, MIRAC) provided the final stop for the two beams of the interferometer. Nulling observations were taken at 10.3, 11.7 and 12.5 $`\mu `$m with 10% bandpass at each wavelength. Images of the science object, HD 100546, were taken at seven different rotations of the interferometer baseline with respect to the sky in order to probe the geometry of the circumstellar dust (discussed in Section 3.1). Thirty-seven sets of 500 frames were taken (2 sets at each combination of wavelength and rotation for all but one of the combinations), each frame with an integration time of 0.5 s. These observations were interlaced with observations of calibrator stars, with the same integration times. Frames were sky subtracted using off-source sky frames taken in between observations. A custom IDL program was used to perform aperture photometry on each image and pick out, for each set of 500 frames, the image with the best instrumental null.

Direct images at 11.7, 18.0 and 24.5 $`\mu `$m were obtained in August 2001 and March 2003. Images were taken with a 3 Hz chop of $`8\mathrm{}`$ in the horizontal direction of the detector, and a nod of $`8\mathrm{}`$ in the vertical direction after every 15 s of integration. The 11.7 $`\mu `$m images were taken to verify extended emission detected by nulling, with a total integration time of 60 s. The 18.0 and 24.5 $`\mu `$m images were taken with the purpose of detecting resolved material at longer wavelengths and characterizing emission from cooler dust. Total integration times were 170 s at 18.0 $`\mu `$m and 210 s at 24.5 $`\mu `$m. Aperture photometry was performed on the direct images of HD 100546 and calibrator stars and relative fluxes for the science object were transformed to absolute fluxes using calibrator fluxes taken from Gezari et al. (1993). The absolute fluxes are 67 Jy at 11.7 $`\mu `$m, 123 Jy at 18.0 $`\mu `$m, and 165 Jy at 24.5 $`\mu `$m, and are good to 10%.

## 3 Results

### 3.1 Nulling Observations

Table 1 shows the source nulls and errors achieved at each combination of wavelength and rotation. The source null for the science object is calculated by subtracting the null achieved on the calibrator star from the instrumental null achieved on the science object, and represents the percentage of light remaining in the image when nulled, compared to the full flux when the pupils are constructively interfered. The best source null is calculated for each set of 500 frames. The adopted values for the source null presented in Table 1 are taken to be the average of the best nulls in each of the two sets of frames taken at each combination of wavelength and rotation. The error in the null is taken to be the difference in the the best nulls in the two sets of frames.

The non-zero source nulls on HD 100546 are indicative of extended emission surrounding the star. Our results show that the source null of the object varies as a function of rotation of interferometer baseline with respect to the sky (Figure A Resolved Circumstellar Disk around the Herbig Ae Star HD 100546 in the Thermal Infrared). At all three wavelengths, the null varies by a factor of 2 or more and appears to have roughly the same dependence. From this dependence, we can infer the presence of an inclined structure, as well as its size and orientation in the following manner. The transmitted signature of the nulling interferometer is an interference pattern with interference fringes along the baseline. If these fringes are parallel to the major axis of a disk, more of the disk’s light will be nulled, resulting in a lower percentage of remaining light. When the fringes are aligned orthogonally to the major axis, the value of the source null will be higher (i.e., there is more light remaining). Therefore, there should be a variation in null with respect to rotation of the baseline. This relation will have a sinusoidal form, with a period of $`180\mathrm{°}`$ and an amplitude dependent on the projected emitting surface of the disk, and consequently its inclination. We have performed a least squares fit of such a function for the expected null (see Fig. A Resolved Circumstellar Disk around the Herbig Ae Star HD 100546 in the Thermal Infrared), $`N=a+bsin(PA+\theta )`$. The value of $`a`$ is a vertical offset (physically related to the size of the disk), $`b`$ is the amplitude of the sine function (related to the inclination of the disk), PA is the position angle of the semimajor axis, and $`\theta `$ is the rotation of the interferometer baseline relative to the sky. We find the parameters shown in Table 2 as the best fits for the data. The values derived for all three wavelengths are in agreement with one another, with average derived uncertainties in the nulls of 6.3, 5.5, and 5.5% for 10.3, 11.7, and 12.5 $`\mu `$m, respectively.

In order to interpret these fits physically, we have calculated the physical parameters which would yield these fit parameters for two simple disk distributions: an inclined Gaussian disk and a ring. Table 3 shows the best fit physical parameters for the nulling data. The best fit position angle of the semimajor axis of the disk is $`140\mathrm{°}`$ to $`160\mathrm{°}`$ (E of N), which is in agreement with the values derived from near-IR coronagraphic studies (Grady et al., 2001; Augereau et al., 2001; Pantin, Waelkens, & Lagage, 2000). The inclination of the disk at 10.3 and 11.7 $`\mu `$m is derived to be $`30\mathrm{°}`$ to $`40\mathrm{°}`$ from face-on, also in agreement with the aforementioned studies, but these fits are only marginally better than a face-on disk when comparing the reduced $`\chi ^2`$ ( = 1.5 and 2.2 for the 10.3 and 11.7 $`\mu `$m fits, respectively) of the fits. The 12.5 $`\mu `$m data show a larger amplitude in the variation of the null, hence the inclination derived is greater, $`60\mathrm{°}`$ from face-on. In this case, there is significant inclination of the disk, as the fit yields reduced $`\chi ^2=2.4`$ as opposed to 7.5 for a face-on disk.

The multiwavelength nature of our observations also allow us to probe differences between the distribution of different species (silicates, PAH, etc.) and the thermal continuum. Our observations at 10.3 $`\mu `$m and 11.7 $`\mu `$m probe emission from silicates and PAH species, while the contiuum emission is roughly probed by the 12.5 $`\mu `$m band (Malfait et al., 1998), although the bandpass of this filter may result in significant emission from PAH and silicates. Our results indicate that the emitting structure is more inclined at 12.5 $`\mu `$m than at the other two wavelengths. This suggests that emission from the thermal continuum may have a more inclined structure than the flux from emission lines of silicates and PAH.

### 3.2 Imaging

The 11.7 $`\mu `$m images verify the presence of resolved emission detected in the nulling data. The images of HD 100546 show an average full-width at half-maximum (FWHM) of $`0.5\mathrm{}`$, while the calibrator star shows a FWHM about 20% smaller. This implies a disk size of about 30 AU, which confirms the disk sizes derived from the nulling data.

The 18.0 and 24.5 $`\mu m`$ direct images show evidence for extended emission as well, with the FWHM values for HD 100546 images on average about 8 - 10% larger than those of the calibrator stars. In order to determine the spatial extent of the extended emission, we constructed an artificial source by convolving an artificial face-on disk signature in the form of a two-dimensional Gaussian, with the PSF from the calibrator star. The artificial image was subtracted from the actual image of HD 100546. The width of the artificial Gaussian disk was varied in steps of 0.1 pix (equivalently $`0.012\mathrm{}`$, or 1.2 AU at 100 pc). We adopt the disk size which resulted in the smallest residual when subtracted from the image of HD 100546. Figure A Resolved Circumstellar Disk around the Herbig Ae Star HD 100546 in the Thermal Infrared shows a typical image before and after subtractions of the artificial source, plotted with the same greyscale. The top-left image shows HD 100546 at 24.5 $`\mu `$m in an unsubtracted image. The top-right image shows the subtracted frame with the smallest residual. The center frames shows the subtraction residuals where the Gaussian disk was about 0.5 pix (FWHM) too small (left), and too large (right). Table 4 shows the results of the model fitting at each wavelength. The disk size adopted for each wavelength is the average of the sizes determined from the observations at the two epochs, with the error bars adopted as the difference in sizes derived for the two epochs.

At these wavelengths, we are probing both the thermal continuum and emission from silicates. As expected from cooler dust, the 24.5 $`\mu `$m emission extends farther out than the 18.0 $`\mu `$m dust. We do note the 24.5 $`\mu `$m band contains a strong emission line from silicates (Malfait et al., 1998) that may contribute significantly on top of the thermal emission. The images also show that there may be evidence for an inclined disk in the subtracted images, as the best subtracted image still shows a roughly symmetric oversubtraction above and below the center of the star. In images at both wavelengths and both epochs, the residuals showed this type of symmetric structure, with peak of the positive residuals on a line orthogonal to the trough of the negative residuals (see center panels, Fig. A Resolved Circumstellar Disk around the Herbig Ae Star HD 100546 in the Thermal Infrared). This structure in the residuals would be expected if the image of HD 100546 was slightly elliptical, perhaps as a result of a resolved inclined disk. We attempted subtractions with artificial sources incorporating an inclined (and rotated) Gaussian disk, rather than a face-on disk. The bottom panel of Figure A Resolved Circumstellar Disk around the Herbig Ae Star HD 100546 in the Thermal Infrared shows the outcome of this subtraction, resulting in a smaller residual than the best subtraction with a face-on disk (top-right panel). The position angle of the artificial disk was changed in steps of $`10\mathrm{°}`$, and the inclination was varied in steps of $`5\mathrm{°}`$ in order to determine the orientation of the artificial disk which resulted in the smallest residual when subtracted from the image of HD 100546. The position angle of the semimajor axis was found to be between $`130\mathrm{°}`$ and $`170\mathrm{°}`$, and the inclination was $`30\mathrm{°}`$ to $`40\mathrm{°}`$ from face-on, all roughly consistent with the orientation derived by previous studies in the near-IR (Augereau et al., 2001; Pantin, Waelkens, & Lagage, 2000; Grady et al., 2001).

## 4 Discussion

We are confident that we have resolved circumstellar emission from HD 100546 at all wavelengths probed in our observations. We wish to compare the physical parameters we have derived for this emission to current models for the circumstellar environments of Herbig Ae stars. Recent models we consider include those described in the Introduction. While observations at all three wavelengths show evidence for an inclined disk, our observations at 10.3 and 11.7 $`\mu `$m are also consistent with a face-on or spherical emitting body. However, the 12.5 $`\mu `$m null variation does provide convincing evidence for an inclined disk. Furthermore, the derived sizes at these wavelengths are increasing with increasing wavelength (equivalently, decreasing temperature), as one might expect. Somewhat puzzling are the derived sizes of the 18.0 and 24.5 $`\mu `$m disks. One would expect the thermal emission at these longer wavelengths to be spatially several times larger than the emission at the shorter wavelengths. If the source of the emission is a continuous flared disk, this relation would be given by T $`r^{0.5}`$ (Chiang & Goldreich, 1997). However, the disks at 18.0 and 24.5 $`\mu `$m are only marginally larger than at shorter wavelengths. This discrepancy suggests that a continuous disk that extends all the way into the dust sublimation radius may not be an accurate model for the dectected emission. Instead we prefer a model with a large inner disk gap, possibly cleared out by a giant protoplanet, as suggested in Bouwman et al. (2003). This would result in the shorter wavelength emission being detected further from the star than expected from a continuous disk, and make the relative sizes of the 10 and 20 $`\mu `$m disks more similar than expected from a T $`r^{0.5}`$ relation.

The dectection of a disk around HD 100546 is also interesting in the context of the observations of Hinz, Hoffmann, & Hora (2001). The previous study performed nulling interferometric observations of three other HAE stars, HD 150193, HD 163296, and HD 179218 and found that none of them had resolved emission. This placed an upper limit on the size of the 10.3 $`\mu `$m disk of 20 AU. This suggests that HD 100546 differs from these other PMS objects, as it appears to have a larger disk at 10.3 $`\mu `$m. A disk such as the one observed around HD 100546 would have been resolvable around the three stars observed by Hinz, Hoffmann, & Hora (2001), and leads us to conclude that the physical structure and/or composition of the circumstellar environment is different in HD 100546 than the HAE stars. This result is consistent with the finding of Bouwman et al. (2003) that the SED of HD 100546 is dissimilar to that of other HAE stars. We also note that the hypothesis of giant protoplanet in the HD 100546 system would provide a natural explanation for the difference between HD 100546 and the stars observed in Hinz, Hoffmann, & Hora (2001). A full analysis of a larger sample of HAE stars is needed to confirm this conclusion, and will be presented in a future paper.

## 5 Acknowledgments

WL acknowledges support from the Michelson Graduate Fellowship. MM and EM acknowledge support through NASA contract 1224768 administered through JPL. EM also thanks the NASA Graduate Student Researchers Program (NGT5-50400) for support. The authors thank the staff of Las Campanas Observatory and the Magellan Project for wonderful support. BLINC was developed under a NASA/JPL grant for TPF. MIRAC is supported by grant AST 96-18850 from the NSF with additional support from SAO.


text:
# Numerical stress response functions of static granular layers

## 1 The simulations

In figure 2, we show a part of a typical simulated layer, where we can see the initial force chain network and the response produced by a localized overload applied to a single grain after the initial forces have been subtracted. The simulations are performed using a classical molecular dynamics (MD) algorithm in three successive stages: preparation, deposition and overloading that are described in the following.

We start with the preparation of $`N`$ polydisperse grains, with radii homogeneously distributed between $`R_{min}`$ and $`R_{max}=2R_{min}`$. These $`N`$ particles are put on the nodes of a grid with aspect ratio $`1:4`$ – this aspect ratio is needed for the appropriate study of the response functions, see below. This lattice ends at its bottom with a fixed horizontal line of similar particles which will be used as the support for deposition – the distances between these bottom grains are small enough to avoid grain evasion.

The deposition stage consists of untying the particles from the grid and, under gravity and horizontal periodic boundary conditions, letting them evolve with random initial velocities. In this ‘rain-like’ deposition the grains interact with elastic and friction forces, and the system evolves under classical Newton equations and molecular dynamics rules (e.g. predictor-corrector, Verlet algorithm). The time to reach the equilibrium depends on the characteristics of the particles. We chose $`\mu =0.5`$ for the contact friction coefficient. The rheology of the contacts is that of Kelvin-Voigt with $`k_n=1000`$ N/m and $`k_t=750`$ N/m for the normal and horizontal contact stiffness. The viscosity $`g_n`$ is chosen in order to get a critical damping. At last, the gravity is set to unity. The equilibrium criteria consist of the five following tests which are applied after each period of 100 MD time steps: (1) the number of gained/lost contacts during this period has to be zero; (2) the number of sliding contacts between particles also has to vanish; (3) the integrated force measured at the bottom of the layer must be equal to the sum of the weight of all the grains; (4) all the particles have to have, at least, two contacts; and (5) the total kinetic energy has to be lower than some low threshold. Once these criteria are all satisfied, the deposition is stopped, and the overloading phase can begin.

The overload is applied to a single grain with a force $`F_0`$ that we wish to be small enough to not cause any rearrangement of the layer structure. The determination of the amplitude of the loading force is then a crucial part of the simulation. For that purpose we made several tests which check the reversibility, the additivity and the linearity of the stress response. These tests (see below for more details) let to conclude that the optimal magnitude of $`F_0`$ is of the order of few times the weight of an average grain. In the overloading procedure, $`F_0`$ is split in 10000 MD time steps, i.e. in each time step the force is increased in 1/10000 of $`F_0`$, and it remains constant and equal to $`F_0`$ after that. Again, the equilibrium criteria are those explicited above, except that the integrated force at the bottom of the layer has to be equal to the weight of all the grains *plus* $`F_0`$.

Only the linearity test is illustrated here, in figure 3: it consists in comparing the force response after loading with different values of $`F_0`$. We find that the response is satisfactorily linear in $`F_0`$ as long as its magnitude does not exceed $`40`$ times the average grain weight. With the other tests we have checked that loading and unloading with the very same $`F_0`$ gives the original contact force distribution with a good precision, and that two simultaneous overloads located at two different positions give a response which is perfectly equivalent to the sum of the two responses computed from the two corresponding single overloads.

## 2 Results and discussion

In order to measure the vertical stress $`\sigma _{zz}`$, we integrate the contact forces over a set of grains. In fact, we do not make use of the usual formula $`\sigma _{\alpha \beta }=1/S_if_\alpha ^ir_\beta ^i`$, where the sum is computed over the contacts between the grains in the ‘volume’ $`S`$ (here we are two-dimensional), $`f_\alpha `$ being $`\alpha ^{th}`$ component of the considered contact force, and $`r_\beta `$ the $`\beta ^{th}`$ component of the corresponding distance vector between the grains in contact. Rather, we take a horizontal ‘line’ of grains of length $`\mathrm{}`$. The corresponding stress $`\sigma _{zz}`$ is then equal to the sum of the vertical components of the forces carried by the contacts of one of the sides of the line – say the upper one – divided by $`\mathrm{}`$. This stress measure can be done at any depth $`z`$ ($`z=h`$ means that the measure is done on the bottom) and centered at any horizontal distance $`x`$ from the overload point. As we are particulary interested in $`\sigma _{zz}`$ profiles along $`x`$ at a given $`z`$, our choice is here better adapted than the usual stress formula which would mix together grains of (slightly) different depths in $`S`$. This is particulary important at small $`z`$. Besides, with this stress definition, the integral of $`\sigma _{zz}`$ over $`x`$ at a fixed $`z`$ is exactly equal to $`F_0`$ at any scale. This property is crucial to normalize and compare data from layers of different thickness.

The response function at depth $`z`$ is obtained by making the difference of the stress profile $`\sigma _{zz}(x,z)`$ measured on the layer with the additional force $`F_0`$, with a copy of the very same layer without overload. This difference is ensemble averaged over several overload positions and layer samples. We present here the results for five different layer sizes, with $`N=100,400,1600,3600`$ and $`6400`$ grains. The $`1:4`$ aspect ratio gives layers of average thickness of $`4,7,14,21`$ and $`29`$ mean particle diameters $`D`$, respectively. The $`10,5,5,4`$ and $`4`$ different realizations give in the end $`62,60,65,84`$ and $`86`$ different overload points – on average, there are six grains between two consecutive overload points.

The results are shown on figure 4. One can see that the normalized profiles measured at the bottom of layers of different thickness $`h`$ collapse when they are linearly rescaled by $`h`$, and show a single peak whose width is of the order of $`h`$. This is precisely the reason why it is important to have layers four times wider than thick – $`2h`$ on both sides of the overload point. When the response is computed in the bulk, the structure of the response profile looks more complicated. The present curve needs however more ensemble averaging and should be taken as preliminary. These results were all obtained using a measure scale $`\mathrm{}=0.77z`$ for all layer sizes. A systematic study of the importance of this integration measure length $`\mathrm{}`$ is under way. We are also working on the calculation of the other stress components $`\sigma _{xx}`$ and $`\sigma _{xz}`$.

These scaling studies are currently extended to other preparation histories, and in particular to the cases of sheared or more anisotropic layers of grains. At last, we plan to test these data against the predictions of anisotropic elasticity which, depending on the values of the different parameters, can give various stress response profiles .

## Acknowledgements

We are indebted to G. Combe for a decise help on the MD codes. We also thank the granular group of the LMDH for stimulating discussions.