text:
# Introduction

## Introduction

The advent of large scale genome projects is transforming the field of molecular evolution (Koonin et al., 2000). The molecular mechanisms of evolution are becoming increasingly amenable to direct observation (Henikoff et al., 1997; Gerstein, 1998: Thornton et al., 1999). and it has become possible to study molecular evolution not only in the context of population genetics, but also by considering the thermodynamic stability of the biomolecules involved in evolution. This “structural” approach has been pioneered by Schuster and co-workers, with a series of studies of neutral networks of RNA secondary structures (Schuster et al., 1994; Huynen et al., 1996; Fontana & Schuster, 1998) and it has been applied to proteins by several groups (Shakhnovich et al., 1996; Bornberg-Bauer, 1997; Bornberg-Bauer & Chan, 1999; Babajide et al., 1997; Govindarajan & Goldstein, 1997, 1998; Bussemaker et al., 1997; Tiana et al., 1998; Mirny & Shakhnovich, 1999; Bastolla et al., 1999, 2000b; Dokholyan & Shakhnovich, 2001). Despite these stimulating studies, however, the structural approach has not yet been used to investigate the classical issues in molecular evolution, as we set out to do here.

In this work we apply to seven of the most studied protein folds the structurally constrained neutral model (SCN) which three of us recently introduced in the context of lattice models (Bastolla et al., 1999; Bastolla et al., 2000b). We compare qualitatively the substitution process obtained from the SCN model and that observed in protein sequence evolution. The SCN model is based on the observation that evolution conserves protein structure much more than protein sequence (Holm & Sander, 1996; Rost, 1997). We assume that all mutations conserving the structure have the same probability of being fixed, thus resulting in a neutral model. The main reason to introduce structure conservation as a working hypothesis is the experimental observation that many mutations do not modify significantly the activity of a protein and its thermodynamic stability, while mutations substantially improving protein functionality are rare (Orencia et al., 2001).

The neutral theory of molecular evolution was introduced in the late 60s by Kimura (1968) and by King and Jukes (1969) to explain the high substitution rates observed in vertebrates as well as the large amount of intra-specific genetic variation observed in most species. According to neutralism, most amino acid substitutions are fixed in the genome of a species not because they bring a selective advantage but due to random genetic drift acting on alleles of equivalent selective value. Kimura’s neutral model predicts that the rate of amino acid substitution of a given protein is approximately constant for different species within major evolutionary groups, independent of the number of individuals and ecology of the species. This was in agreement with earlier observations on protein evolution which lead to postulate a kind of “molecular clock” (Zuckerkandl & Pauling, 1962).

According to Kimura’s theory, the fraction of amino acids neutrally substituted in a time $`T`$ is a Poissonian variable of expectation value $`kT`$ where $`k`$, the substitution rate, is different for different proteins but does not change for different species. A subsequent study by Ohta and Kimura (1971) measured the variance of the substitution process acting on different species, finding that it is larger than the mean, i.e. the process is non-Poissonian. Such a result was confirmed by the more sophisticated analyses by Langley and Fitch (1973) and Gillespie (1989). This and other observations lead first Ohta and then Kimura to adopt, in place of the original neutral model, a model based on slightly deleterious mutations (Ohta, 1976), and Gillespie to reject in toto the neutral theory favouring the hypothesis that most substitutions in protein sequences are fixed by positive selection (Gillespie, 1991). Takahata, however, showed that an extension of the neutral theory, the fluctuating neutral space model (Takahata, 1987), accounts for the non-Poissonian statistics of substitutions.

One of the goals of this study is to investigate the consequences of structural conservation on the properties of neutral networks and on the substitution process associated to them. We show that neutral evolution does not lead to a Poissonian substitution process. This result complements the fluctuating neutral space model by Takahata (1987), and suggests that arguments against the neutral theory based on the fact that the substitution statistics is non-Poissonian (Gillespie, 1991) are not conclusive. A deeper understanding of the mechanism of neutral evolution will help to single out the perhaps less common but more interesting cases of positive selection as, for instance, functional changes and responses to changes in the environment. It can also be useful for calibrating the molecular clocks used to reconstruct phylogenetic trees, whose reliability is severely limited by the fluctuations of the substitution rate (Ayala, 1997).

Another interesting application of the SCN model is the possibility to distinguish between functional and structural conservation. By simulating neutral evolution we identify the key positions which are more difficult to mutate. We identify them as structurally conserved positions, and those positions conserved in actual evolution but not in the SCN model as functionally conserved ones. Practically all residues whose functional role is known belong to this class. Most positions are not conserved in the SCN model, as similarly observed in actual evolution data. We identify them as neutrally evolving positions, and argue that their preeminence is an evidence of the importance of neutral evolution. Finally, a small number of positions appear structurally important in the SCN model but are not significantly conserved. This could be due to a limitation either of the SCN model or of the protein database, but it could also be a clue of structural changes, possibly positively selected. Other methods to identify computationally structurally important positions have been proposed recently (Kannan & Vishveshwara, 1999; Cecconi et al., 2001). In particular, another method based on simulated evolution has appeared in a recent preprint after this work had been completed (Dokholyan & Shakhnovich, 2001).

As most computational studies of protein evolution, the SCN model is based on an approximate stability criterion relying on the $`Z`$-score (Bowie et al., 1991; Goldstein et al., 1992) and on a folding parameter measuring the degree of correlation of the energy landscape (Bastolla et al., 1999). While these parameters can not predict precisely the thermodynamic stability of a specific sequence, our previous studies show that they correlate with the observed stability. Thus, we expect that the statistical properties derived from the analysis of a large number of sequences capture real features of protein evolution.

## Structurally Constrained Neutral Model

Following Kimura, we divide mutations in two classes: those which result in inactivating the protein, which are regarded as lethal and can not spread in the population, and those after which the protein remains active, which are regarded as selectively neutral. This mutational spectrum implies that protein sequences evolve on a neutral network, i.e. a set of sequences where the protein is active and which can be connected through point mutations. Under this mutational spectrum, fixation of slightly deleterious mutations can not take place, since these are not included in the model, as well as advantageous mutations. This is of course an important limitation of neutral models.

Kimura’s neutral model assumes that the rate of appearence of neutral mutations is constant throughout evolution. In a paper of 1977, however, Kimura comments that rate constancy may not hold exactly (Kimura, 1977). Our model is not based on any assumption on the neutral mutation rate. Rather, we compute the effect of mutations on protein stability using an effective model of protein folding (Bastolla et al., 2000a) which provides us with a genotype to phenotype mapping. In this respect, the rate of occurrence of neutral mutations is an outcome of the model. It turns out that this rate shows very broad fluctuations throughout evolution. As we shall see below, the variance in evolutionary rates predicted by our model is in qualitative agreement with observations of protein evolution (Gillespie, 1991).

The model does not take into account population dynamics. This is based on the fact that, within Kimura’s model, the substitution rate is not influenced by population size. An extension of Kimura’s model to take into account small variations in the neutral mutation rate confirmed this result (Bastolla & Peliti, 1991). However, population size might influence the substitution rate if the rate of neutral mutations shows broad fluctuations, as observed here. The explicit inclusion of population genetics into the model would be needed to investigate this interesting possibility.

A neutral network is defined starting from a protein sequence in the Protein Data Bank (PDB). The corresponding protein structure has to remain thermodynamically stable during evolution. Thermodynamic stability is evaluated through an effective model of protein folding. The folding parameters of the native structure, computed through the model (see Materials and Methods), must be above 98.5 percent of the value they have in the PDB sequence. Sequences where this condition is met are named viable sequences. A neutral network is a set of viable sequences which can be connected to the starting sequence through point mutations passing on other viable sequences. Thus sequences on a neutral network share the same protein fold and are evolutionarily connected. For every amino acid sequence $`𝐀`$ in the neutral network we can measure the fraction of neutral neighbors $`x(𝐀)`$, which is the fraction of its possible point mutations which are viable.

We model protein evolution at the level of a single sequence. During evolution, the sequence moves on the neutral network generating an evolutionary trajectory, i.e. a list of subsequently visited sequences belonging to the neutral network. In the present context, the only relevant quantity is their fraction of neutral neighbors $`x(𝐀)(0,1]`$. Thus an evolutionary trajectory can be represented through a very long list $`𝐱=\{x_1,x_2,\mathrm{}\}`$.

The rate of occurrence of a neutral mutation starting from sequence $`𝐀`$ is the product of a constant mutation rate times the probability $`x(𝐀)`$ that the mutation is viable. Thus the neutral mutation rate is not constant in the framework of our model, and can be explicitly computed by computing $`x(𝐀)`$ for all sequences in an evolutionary trajectory. The statistics of the substitution process can then be obtained by coupling the evolutionary trajectory generated as above to a Poissonian mutation process according to the following rules:

1. Mutation process: The number of mutations in a time $`t`$ is a Poissonian variable of average $`\mu t`$.
2. Acceptance process: Given one realization of the evolutionary trajectory and a number of mutations $`k`$, the conditional probability that $`n`$ of them are accepted is the product of $`n+1`$ geometric distributions of parameters $`1x_i`$:

$$P_{\mathrm{acc}}\left(nk\right)=\left(\underset{i=1}{\overset{n}{}}x_i\right)\underset{\{m_j\}}{}\underset{j=1}{\overset{n+1}{}}(1x_j)^{m_j},$$
(1)
where the $`\{m_j\}`$ are all integer numbers between zero and $`kn`$ satisfying $`_{j=1}^{n+1}m_j=kn`$. In other words, the probability that a mutation is accepted is $`x_1=x(𝐀_1)`$ as long as the protein sequence is $`𝐀_1`$, $`x_2=x(𝐀_2)`$ as long as the sequence is $`𝐀_2`$, and so on.

Two kinds of random variables must be distinguished. We indicate by angular brackets the average over mutation and acceptance process for a given realization of the evolutionary trajectory and by an overline the average over evolutionary trajectories. The variance of the substitution process can be decomposed in two components:

$`V(S_t)`$ $`=`$ $`V_\mu (S_t)+V_x(S_t)`$ (2)
$`=`$ $`\left(\overline{S_t^2}\overline{S_t^2}\right)+\left(\overline{S_t^2}\overline{S_t}^2\right).`$ (3)
The first term, $`V_\mu `$, is the variance of the mutation and acceptance process, averaged over evolutionary trajectories. The second term, $`V_x`$, is the variance of the substitution rate with respect to different evolutionary trajectories. This term, which is not present in the standard neutral model, explains why the variance of the number of substitution is typically larger than its mean value, contrasting with a Poissonian process.

If all sequences have the same fraction of neutral neighbors $`x(𝐀)x`$, the number of substitutions in a branch of length $`T`$ is Poissonian with mean $`\mu Tx`$ and the substitution rate is equal to $`\mu x`$ as in Kimura’s model. If $`V_x`$ is not zero, the substitution distribution is more complicated and has to be computed numerically using the evolutionary trajectories simulated.

## Numerical results

### Folding of random sequences

As a preliminary analysis, we measured the distribution of the folding parameters $`\alpha `$ and $`Z^{}`$ (see Materials and methods) of the native structures considered in this work for random sequences of the same length of the corresponding PDB sequences. On over 20,000 attempts, we always found folding parameters much lower than for PDB sequences. The only exception was the smallest protein, the 53-residues rubredoxin, for which a single random sequence had one of the stability parameters comparable to that of the PDB sequence, even if still smaller than it. This result is consistent with the work of Keefe and Szostak who were recently able to select ATP binding proteins from a random library of more than $`10^4`$ sequences (Keefe & Szostak, 2001). We note that it is possible to evaluate the size of the neutral network from the joint distribution of $`\alpha `$ and $`Z^{}`$ , but to this end better statistics are needed than those obtained here.

### Connectivity of neutral networks

For sequences $`𝐀`$ belonging to the neutral network, the fraction of neutral neighbors $`x(𝐀)`$ counts the fraction of all possible point mutations of $`𝐀`$ which still fall into the neutral network. We measured this quantity for at least 20,000 sequences for each fold, finding that it has a broad distribution (see Fig. 1). The shape of the distribution is qualitatively similar for all of the studied proteins, but in the case of cytochrome c, the distribution is shifted to lower connectivities. Results for all seven folds are summarized in Table 1.

The connectivity landscape $`x(𝐀)`$ is locally correlated for distances in sequence space of the order of at least ten substitutions. The correlation function $`C(t)=\left(x(𝐀_{t_0})x(𝐀_{t_0+t})x(𝐀)^2\right)/\sigma _x^2`$ between connectivities of two sequences at distance of $`t`$ steps decays similarly for all proteins and can be fitted to a stretched exponential $`C(t)\mathrm{exp}\left((t/\tau )^\eta \right)`$, with exponents $`\eta `$ ranging from $`0.60`$ to $`0.66`$ and correlation lengths $`\tau `$ ranging from $`1.8`$ to $`2.8`$. Thus correlations decay to one tenth after about ten substitutions (data not shown).

### Substitution process

The broadness of the connectivity distribution directly implies that the substitutions process fluctuates more than a Poissonian process, i.e. it is overdispersed. We computed average and variance of the substitution process numerically, using the evolutionary trajectories generated in our simulations. Results for myoglobin are shown in Fig. 2. Notice that the substitution rate $`\overline{S_t}/t`$ is roughly constant in time, and the total dispersion index $`R(t)=\overline{V(S_t)}/\overline{S_t}`$ takes values between 1.0 at small $`t`$ and 1.9 at large $`t`$, consistent with the value $`R=1.7`$ estimated by Kimura for myoglobin.

Due to local correlations in sequence space, different evolutionary trajectories $`x(𝐀_1)\mathrm{}x(𝐀_n)`$, representing different populations, give different mean and variance of the substitution process over short time scales. This phenomenon produces new lineage effects, i.e. apparently varying substitution rates in different branches of the phylogenetic tree. To illustrate them, we show in Fig. 3 the mean $`S_t(\{x\})`$ and the variance $`V(S_t,\{x\})`$ for three realizations of the evolutionary trajectory $`\{x\}`$. Such an effect could overshadow the generation time effect for replacement substitutions (Britten, 1986; Li et al., 1987; Gillespie, 1991). It could also been responsible for the wide fluctuations in the substitution rate for different lineages observed by Ayala (1997).

### Sequence similarity

One may ask whether sequences sharing the same fold must have a high level of similarity. To investigate this question, we measured the distribution of sequence similarity for sequences in the neutral network obtained through simulation of our model as well as for homologous sequences in public databases. Similarity between two aligned sequences of the same length is defined as the fraction of positions where the same amino acid appears. Results are showed in Fig. 4 for the globin fold. Similar results have been obtained for all other folds.

Sequences in the neutral network (solid line) have an average similarity only slightly larger than random sequences, for which a Gaussian-like distribution of average value $`1/20`$ is expected. This confirms a previous finding of three of us for neutral networks of lattice structures (Bastolla et al., 1999). As previously observed by Rost (1997), the same result holds also for sequences in the FSSP family of sequences sharing the same structure (Holm & Sander, 1996), which are showed as dotted line. Low similarity for sequences with the same fold has also been found in a recent computational study based on sequence optimization for native protein structures (Dokholyan & Shakhnovich, 2001). The dashed line shows the similarity distribution for sequences in the PFAM database of similar sequences (Bateman et al., 2000). The PFAM family has on the average a much larger similarity. This is in part due to the fact that in this case sequence similarity must be large enough for the homology to be detected and in part to the fact that proteins in a PFAM family are subject to stronger functional conservation than proteins in the FSSP family. For the globin family, which has been intensively studied, even the PFAM similarity is very low.

### Residue conservation

The SCN model identifies structurally conserved residues. In this respect, the present results can be compared to the model of Shakhnovich and coworkers (Shakhnovich et al., 1996; Dokholyan & Shakhnovich, 2001) based on sequence optimization (Shakhnovich & Gutin, 1993; Shakhnovich, 1994), and to the bioinformatic studies of Ptitsyn’s group (Ptitsyn, 1998; Ptitsyn & Ting, 1999).

We evaluate the conservation of each position measuring its rigidity (see Materials and Methods). This was done for each fold using three different sets of sequences: (i) Sequences obtained from simulations of our neutral model; (ii) Homologous sequences in the PFAM database (Bateman et al., 2000); (iii) Sequences with the same structure in the FSSP database (Holm & Sander, 1996). The PFAM family often contains orthologous proteins performing the same function while in the FSSP family different functions may be present and function conservation plays a less important role. Nevertheless, there is usually a good correlation between the rigidity of a given position evaluated through the PFAM and FSSP databases. An exception is the TIM barrel family, one of the most common folds, used to perform different functions, each approximately associated to a different PFAM family. In this case, the two sets of rigidities show no correlation. For this fold, using a bioinformatic analysis, Mirny and Shakhnovich (1999) found evidence of functional conservation (the same functional positions tend to be conserved in all functional families, although with different residues in each functional family) but could not find evidence of structural conservation.

Before turning to the analysis of conserved positions, we observe that there are three reasons why sequence databases may tend to overestimate structurally based conservation. The first one is the small size of databases. The second one is the fact that many sequences are evolutionarily related: databases usually provide biased samples of the tree of life. To reduce these effect, we limit our analysis to sequences that do not have similarity larger than a threshold $`q_{\mathrm{thr}}`$ which we choose equal to 0.85 in order not to reduce too much database size, and try to estimate the maximal conservation that one would observe with a database of similar size and correlations, in the null hypothesis that all positions are equivalent. The third reason is that many residues are conserved on functional grounds, sometimes even in the FSSP database, and it may be difficult to distinguish them from structurally conserved residues.

Conservation in the neutral network only expresses structural conservation, thus the comparison between rigidities predicted by the SCN model and observed in evolution may allow to single out functionally conserved positions or positions involved in interactions with cofactors, which are not represented in our model. We tested this for two well studied protein families: the globin family and the cytochrome c family. In both cases structurally conserved positions identified by the SCN model coincide with those identified in previous bioinformatic studies as part of the folding nucleus (Ptitsyn, 1998; Ptitsyn & Ting, 1999), and additional structurally conserved positions are found. For other protein families less is known about functional residues, but the few ones which are identified in the SwissProt file are recognized as such by the SCN model. In our analysis of the globin family, positions in contact with the heme are not regarded as structural, even if the heme plays also an important stabilizing role, since interactions between amino acids and cofactors are not considered in the model, and they are much more specific than interactions between amino acids.

In Fig. 5 we compare the rigidities obtained from our model to those measured in the FSSP family for the myoglobin fold. Each point represents a position on the native structure. The dotted lines in the figure are rough estimates of the maximal rigidity expected in a random situation, i.e. all equivalent residues and same distribution of similarity as in the set examined. Only residues more rigid than that are considered significantly conserved. Many of the most conserved residues are in contact with the heme group (large circles). A notable exception is Pro37<sup>*</sup><sup>*</sup>* Residues are labeled in the order in which they are listed in the PDB file of the structure 1a6g, which is strongly conserved and not in contact with the Heme group. Although the conservation of this residue has not been fully explained so far, Ptitsyn and Ting report that it may be due to functional reasons (Ptitsyn & Ting, 1999). The three positions most conserved according to the SCN model coincide with structural positions identified in the bioinformatic analysis by Ptitsyn and Ting. They are, in order of rigidity, Leu115, Trp14, Met131. Val10 is rather conserved both in our model and in the bioinformatic study. The remaining two positions identified by Ptitsyn, Ile111 and Leu135, are not among the most conserved in our model, although they are above the average. In addition, there are eight more positions significantly conserved in the SCN, whose evolutionary conservation is somewhat lower (but above the average in all cases except His119). They are: Val13, Val17, His24, Leu69, Leu76, His119, Phe123, Ala134. Interestingly, structurally conserved residues form a cluster, so that it has been proposed that they play the role of “folding nuclei” (Shakhnovich et al., 1996; Ptitsyn 1998; Mirny & Shakhnovich, 2001). A similar situation applies also for the case of cytochrome c: two of the positions identified in (Ptitsyn, 1998) are the most conserved in our model (Phe7, Leu74), another one is significantly conserved (Trp77) and the fourth one is not present in the structure we choose as reference (PDB code 451c). Moreover, there are three positions significantly conserved in our model and in the FSSP alignment (Tyr27, Ile48, Val66) and one conserved in our model but not in the alignment (Gly36).

We show a similar plot also for rubredoxin, a small bacterial protein involved in electron transport. In this case, we studied two homologous proteins, one from a mesophylic and one from a thermophylic bacterium. Their sequences have 57% similarity and belong to the same PFAM and FSSP classes. Although the structures are rather similar, the stability of the thermophylic protein, as measured by the $`Z^{}`$ and $`\alpha `$ parameters, is higher than the stability of the mesophylic protein, as it should be. This result supports our choice of the stability criteria. We compare the rigidities obtained from the SCN model for the two structures in Fig. 7 (upper panel). There is a remarkable correlation, despite the fact that results are obtained from independent evolutionary runs with different selection parameters. In Fig. 7 (lower panel) we compare the rigidity observed in the SCN model (for mesophylic rubredoxin) with the rigidities observed in the PFAM and FSSP databases.

Finally, we show in Fig. 8 a scatter plot of rigidities for the protein showing the worst correlation between predicted and observed rigidities: the TIM barrel, one of the most common folds, used for several enzymatic functions. The one that we studied is a triose phosphate isomerase functioning in the glycolysis. In this case, there is also no observable correlation between rigidities in the PFAM and FSSP databases, and rigidities in the FSSP class are very low, in particular because several residues are deleted in many sequences.

Unlike for other folds, in the case of ubiquitin the structurally conserved positions are distributed along the main axis of the protein. There is also a conserved polar position on one loop. For lysozyme as well the most conserved residues form a non compact core. In all other cases the structurally conserved residues form a hydrophobic cluster which is rather compact. There is some correlation between conservation and burial in the interior of the protein as measured by the number of contacts, but burial alone does not explain all of the conservation. We notice that also in our previous lattice simulation the most conserved positions were those in the interior of the lattice structure (Bastolla et al., 1999).

## Discussion

In this work we studied a neutral model of protein evolution based on structure conservation. For all of the examined folds, local connectivities of neutral networks are broadly distributed. This result implies that even in neutral evolution the number of substitutions fluctuates more than a Poissonian variable, i.e. it is overdispersed. Therefore overdispersion can not by itself be used as a test for positive selection, as argued for instance by Gillespie (1991). Our results, nevertheless, show that the dispersion index of the SCN substitution process is typically small, rarely overcoming four. Thus, proteins with extremely high dispersion index, as some of those studied by Gillespie (1989) or Ayala (1997) are not likely to have evolved in this way. The average substitution rate is almost constant in time, but it may vary considerably for different evolutionary trajectories, corresponding to different branches of the phylogenetic tree. This fact should be taken into account when studying lineage effects such as the generation time effect (Britten, 1986; Li et al., 1987; Gillespie, 1991).

By simulating neutral evolution, we identified structurally conserved positions and compared them to evolutionarily conserved positions in known protein families. The comparison is favorable for myoglobin, cytochrome c, lysozyme, ribonuclease and rubredoxin, while for ubiquitin and the TIM barrel correlation between predicted and observed conservation is almost absent. We note, however, that the TIM barrel shows very little structural conservation, and the small size of the ubiquitin family makes the comparison not conclusive. The plots comparing conservation in simulated evolution (on the abscissa) to conservation in real evolution (on the ordinate) can be divided in four parts. In the upper left quadrant there are positions not conserved in the SCN model but conserved in evolution. We suggest that most of them are conserved for functional reasons or because of interactions with cofactors, which are not taken into account in our protein model. Positions of known functional importance belong to this class, but not enough is known on protein function to prove our interpretation in all cases. In the upper right quadrant there are positions conserved both in the SCN model and in the databases, whose conservation is likely to have a structural ground. For this small subset the rigidities that we predict are correlated to the observed ones. Interestingly, those positions form spatial clusters which have been identified with folding nuclei (Shakhnovich et al., 1996; Ptitsyn, 1998; Mirny & Shakhnovich, 2001). Although we can not discuss such interpretation, since our evolutionary algorithm does not take into account folding kinetics, it is to be expected that positions important for stability also play an important kinetic role. In the bottom left quadrant there are positions not conserved neither in the SCN model nor in the databases. These positions are likely to be the main actors in neutral evolution. Last, in the bottom right quadrant there are few positions conserved according to the SCN model which do not appear to be evolutionarily conserved. Barring artifacts due to the SCN model, we should consider the possibility of conservation with low rigidity (but typically much higher than random). In order to verify whether this is the case, we need larger and less correlated protein classes. Another possibility is that these positions are frequently substituted because they can produce structural changes, possibly positively selected. Although this possibility is rather speculative, it would be interesting to investigate it in more detail.

Our results are based on an approximate stability criterion relying on the $`Z`$-score (Bowie et al., 1991; Goldstein et al., 1992) and on a parameter measuring the degree of correlation of the energy landscape (Bastolla et al., 1999). While such criterion may not be suitable for the quantitative prediction of the thermodynamic stability of particular proteins, we believe that the statistical properties of the SCN model reflect those of actual protein evolution. This confidence has several grounds. First, in the study of a lattice model, three of us have previously applied a rigorous criterion of stability, and compared it to a criterion obtained from the $`Z`$-score (Bastolla et al., 2000b). Although the two criteria give different responses for specific sequences, it is possible to choose a threshold such that most sequences selected with the $`Z`$-score criterion are also selected with the rigorous criterion. Second, the present results are robust with respect to changes in the selection thresholds and stability criteria. Third, we tested our stability parameters on a large number of sequences obtained from mutations of a TIM barrel enzyme, whose phenotypic effect has been experimentally measured in a recent paper (Silverman et al., 2001). We found that, even if our criterion can not predict the effect of individual mutations, the latter is correlated to the $`\alpha `$ parameter with correlation coefficient 0.4.

The present results can be used in the rational approach to directed evolution of biocatalysts (Altamirano et al., 2000; Voigt et al., 2001) since we identify sites that are more tolerant to mutations and therefore likely targets for evolutionary improvement. This is a remarkable possibility, since it indicates how results based on the assumption of neutral evolution can be used to search for positive substitutions.

## Materials and Methods

### Protein model

We represent a protein structure by its contact matrix $`C_{ij}`$, where $`C_{ij}=1`$ if residues $`i`$ and $`j`$ are in contact and $`C_{ij}=0`$ otherwise. Two residues are considered in contact if any two of their heavy atoms are closer than 4.5 Å. The effective free energy associated to a sequence of amino acids $`𝐀`$ in the configuration $`𝐂`$ is approximated as a sum of pairwise contact interactions,

$$E(𝐀,𝐂)=\underset{i<j}{}C_{ij}U(A_i,A_j),$$
(4)
where $`A_i`$ labels one of the twenty amino acid types and $`U(a,b)`$ is a $`20\times 20`$ symmetric interaction matrix. Here we use the matrix derived by Bastolla et al. (2000a), which describes accurately the thermodynamic stability of a large set of monomeric proteins (Bastolla et al., 2001).

Three remarks are needed: First, the effective energy parameters implicitly take into account the effect of the solvent and depend on temperature. They express free energy rather than energy. Second, the effective energy of a structure is defined with respect to a completely extended reference structure where no contacts are formed and which sets the zero of the energy scale. Third, one can derive from the database not the parameters $`U(a,b)`$ themselves but the dimensionless quantities $`U(a,b)/k_BT`$. It is thus important to use dimensionless parameters to evaluate the stability of the protein model.

### Candidate structures

We generate candidate structures for a protein sequence of $`N`$ residues by generating all possible gapless alignments of the sequence with structures in the Protein Data Bank. This procedure is called threading. In this way, we typically generate several hundreds of thousands of protein-like structures per sequence. In the present context, threading is directly used to produce the contact maps of the candidate structures. In order to speed up the computations, we use a non redundant subset of the PDB excluding proteins with homologous sequences, selected by Hobohm & Sander (1994).

### The folding parameter $`\alpha `$

For a given sequence $`𝐀`$, the energy landscape is well correlated if all configurations of low energy are very similar to the configuration of minimal effective energy, $`C^{}(𝐀)`$. Structure similarity is measured by the overlap $`q(𝐂,𝐂^{})`$, counting the number of contacts that two structures have in common and normalizing it through the maximal number of contacts, so that $`q`$ is comprised between zero and one. In a well correlated energy landscape, the inequality holds

$$\frac{E(𝐀,𝐂)E(𝐀,𝐂^{})}{|E(𝐀,𝐂^{})|}\alpha (𝐀)\left(1q(𝐂,𝐂^{})\right),$$
(5)
stating that the energy gap between each alternative structure $`𝐂`$ and the ground state $`𝐂^{}`$, measured in units of the ground state energy, is larger than a quantity $`\alpha (𝐀)`$ times the structural distance $`1q(𝐂,𝐂^{})`$. The dimensionless quantity $`\alpha (𝐀)`$, which is the largest quantity for which the above inequality holds, can be used to evaluate the folding properties of sequence $`𝐀`$. For random sequences, many different configurations have quite similar energy and $`\alpha (𝐀)0`$. In this case the energy landscape is rugged, the folding kinetics is very slow and the thermodynamic stability with respect to variations in the solvent is very low. In contrast, computer simulations of well designed sequences have shown that, when $`\alpha (𝐀)`$ is finite, the folding kinetics is fast and the stability with respect to changes in the energy parameters as well as mutations in the sequence is very high.

Our algorithm computes the parameter $`\alpha (𝐀)`$ for a fixed target configuration $`𝐂^{}`$ and a large number of sequences $`𝐀`$. We thus indicate this parameter as $`\alpha (𝐀,𝐂^{})`$, since we do not know a priori that $`𝐂^{}`$ has lowest energy. Notice however that, if $`\alpha (𝐀,𝐂^{})`$ is positive, all alternative structures have higher energy than $`𝐂^{}`$. We impose that $`\alpha (𝐀,𝐂^{})`$ is larger than a positive threshold $`\alpha _{\mathrm{thr}}`$ for sequences $`𝐀`$ belonging to the neutral network.

### $`Z`$-score

The Z-score $`Z(𝐀,𝐂^{})`$ (Bowie et al., 1991; Goldstein et al., 1992) is a measure of the compatibility between a sequence $`𝐀`$ and a structure $`𝐂^{}`$, widely used in structure prediction. It depends on an effective energy function, and measures the difference between the energy of sequence $`𝐀`$ in configuration $`𝐂^{}`$ and its average energy in a set of alternative configurations, $`\{𝐂\}`$, in units of the standard deviation of the energy:

$$Z(𝐀,𝐂^{})=\frac{E(𝐀,𝐂^{})E(𝐀,𝐂)}{\sqrt{E(𝐀,𝐂)^2E(𝐀,𝐂)^2}}.$$
(6)
When sequence $`𝐀`$ folds in structure $`𝐂^{}`$ their corresponding Z score is very negative.

Given the above definition, one has still to specify how to select alternative structures. A possibility, often used for lattice models (Mirny & Shakhnovich, 1996) is to assume that alternative structures are maximally compact, randomly chosen structures, whose average energy can be estimated as $`E(𝐀,𝐂)_𝐂=Nc_{\mathrm{max}}e(𝐀)`$. Here, $`Nc_{\mathrm{max}}`$ is the maximal number of contacts of candidate structures and $`e(𝐀)`$ is the average energy of a contact, averaged over all possible contacts formed by sequence $`𝐀`$. This leads to introduce the parameter

$$Z^{}=\frac{E(𝐀,𝐂^{})/Nc_{\mathrm{max}}e(𝐀)}{\sqrt{e^2(𝐀)e(𝐀)^2}}.$$
(7)
The use of $`Z^{}`$ has two main advantages: First, it makes the value of the $`Z`$-score much less sensitive to chain length $`N`$ and to the particular set of alternative structures used; second, the evaluation of $`Z^{}`$ is much faster than that of the $`Z`$-score. This is necessary in order to explore efficiently sequence space.

### Sampling the neutral network

Our algorithm explores the neutral network of a given protein starting from its PDB sequence $`𝐀_0`$ and iterating the following procedure: At time step $`t`$, (i) The number of viable neighbors of sequence $`𝐀_t`$ is computed; (ii) The sequence $`𝐀_{t+1}`$ is extracted at random among all the viable neighbors of $`𝐀_t`$. In this way we generate a stochastic process along the neutral network which simulate neutral evolution and looses memory of the initial sequence very fast.

Sequence $`𝐀`$ is regarded as viable if both parameters $`\alpha (𝐀,𝐂^{})`$ and $`Z^{}(𝐀,𝐂^{})`$ are above predetermined thresholds, chosen as 98.5 percent of the values of those parameters for the sequence in the PDB. This enforces conservation of the thermodynamic stability and folding capability of the native structure $`𝐂^{}`$. We verified that the observed behavior does not change qualitatively for thresholds between 95% and about 100% of the PDB values.

We impose strict conservation of the cysteine residues in the PDB sequence, and do not allow any residue to mutate to cysteine, since a mutation changing the number of cysteine residues would leave the protein with a very reactive impaired cysteine that would probably affect its functionality. Accordingly, the total number of neighbors tested is $`X_{\mathrm{tot}}=18(NN_{\mathrm{cys}})`$, where $`N`$ is the number of residues and $`N_{\mathrm{cys}}`$ is the number of cysteine residues in the starting sequence.

The total number of viable point mutations, $`X(𝐀)`$, expresses the local connectivity of the neutral network. We normalize it by the total number of neighbors, $`X_{\mathrm{tot}}`$, getting the fraction of neutral neighbors, $`x(𝐀)=X(𝐀)/X_{\mathrm{tot}}(0,1]`$.

To compute $`x(𝐀)`$, we have to evaluate the $`\alpha `$ parameter for all sequences $`𝐀^{}`$ obtained through a point mutation of sequence $`𝐀`$. From Eq.(5) we note that the $`\alpha `$ parameter can be obtained from the configuration with the highest destabilizing power, i.e. the highest value of the energy gap divided by the structural distance from the native configuration. These change from sequence to sequence, but it is expected not to change very much for neighboring sequences. Thus, in order to speed up the computation of $`\alpha (𝐀^{})`$, instead of considering all candidate configurations we consider only the 50 configurations with the highest destabilizing power (i.e. the energy gap divided by the structural distance from the native configuration) for sequence $`𝐀`$ and compute their mutated destabilizing power using sequence $`𝐀^{}`$. The $`\alpha `$ parameter is then obtained from the configuration with the highest destabilizing power. This procedure could slightly overestimate $`\alpha (𝐀^{})`$ since not all configurations are used, but we have checked that the error introduced in the $`x`$ value is in all cases below 0.1%.

### Substitution process

Given an evolutionary trajectory $`\{x_1,x_2,\mathrm{}\}`$, the distribution of the number of substitutions taking place in a time $`T`$ can be computed by considering Eq.(1), where $`k`$, the number of attempted mutations, is a Poissonian variable of average value $`\mu T`$.

In order to handle the computation, we divide all values of $`x_i`$ in $`M`$ classes, choosing $`X_a`$ as representative value of all $`x_i`$’s belonging to class $`a`$. The number of operations needed to evaluate the substitution probability increases exponentially with the number of classes $`M`$. At the same time, the evaluation becomes more and more accurate as $`M`$ increases. We chose $`M=6`$ in our numerical computations as a reasonable compromise between accuracy and rapidity, checking that larger values of $`M`$ introduce only small changes.

### Rigidity

A measure of the conservation profile for a set of evolutionarily related sequences can be obtained measuring the rigidity of each position $`i`$,

$$R(i)=\underset{a}{}f_i(a)^2,$$
(8)
where $`f_i(a)`$, $`a=1,\mathrm{}20`$ is the frequency with which amino acid $`a`$ is observed at position $`i`$, normalized so that $`_af_i(a)1`$. Deletion of position $`i`$ in a sequence is regarded formally as a 21st amino acid. A large rigidity $`R(i)`$ means that position $`i`$ is highly conserved. For unconstrained positions and in absence of deletions, $`f_i(a)=1/M`$, where $`M`$ is the number of amino acids, and $`R(i)=1/M`$. In general, rigidities are larger than $`1/M`$ because of the finite size of the sequence set and because sequences in the set are correlated due to common evolutionary origin. Since cysteine residues are strictly conserved, we always get $`R(i)=1`$ for them. Thus we omit these residues from the analysis of conservation.

### PFAM and FSSP databases

We compare the rigidity measured in the set of neutral sequences generated with the present method with the rigidity obtained from two databases: the PFAM database (Bateman et al., 2000) and the FSSP database (Holm & Sander, 1996). The PFAM database is a collection of families of homologous sequences obtained by multiple alignment. Since multiple alignment methods work only for sufficiently high similarity, there are no sequences of low similarity in this database. The FSSP database is a collection of protein classes sharing the same fold (as determined by the program of structural alignment DALI (Holm & Sander, 1996)). Since the structures must be experimentally known, the FSSP database is usually smaller than the PFAM database. However it includes in the same class distant homologs whose evolutionary relationship can not be detected by means of sequence comparison alone. Due to database biases, many sequences in the PFAM and FSSP databases are highly similar. To reduce this effect, proteins with similarity higher than a threshold $`q_{\mathrm{thr}}=0.85`$ to any other protein have been eliminated.

## Acknowledgments

UB thanks for interesting discussions Bill Eaton, Carlos Briones, Nikolai Dokholyan, Leo Mirny, Victor Parro, Burkhard Rost, Alfonso Valencia and Roeland Van Ham, and the Institut of Crystallography of the Free University of Berlin for logistic help. UB is supported by INTA (Spain). MP is supported by the Max Planck Gesellschaft (Germany). HER is supported by INFN (Italy). MV is supported by EMBO and by the Royal Society (UK).

Abbreviations: PDB Protein Data Bank, SNC Structurally Constrained Neutral Model, FSSP Fold classification based on Structure-Structure alignment of Proteins.


format: text

text:
# Monte Carlo transition probabilities

## 1 Introduction

When Monte Carlo methods are used to compute the spectra of astronomical sources, it is advantageous to work with indivisible monochromatic packets of radiant energy and to impose the constraint that, when interacting with matter, their energy is conserved in the co-moving frame. The first of these constraints leads to simple code and the second facilitates convergence to an accurate temperature stratification.

For a static atmosphere, the energy-conservation constraint automatically gives a divergence-free radiative flux even when the temperature stratification differs from the radiative equilibrium solution. A remarkable consequence is that the simple $`\mathrm{\Lambda }`$-iteration device of adjusting the temperature to bring the matter into thermal equilibrium with the Monte Carlo radiation field results in rapid convergence to the close neighbourhood of the radiative equilibrium solution (Lucy 1999a). An especially notable aspect of this success is that this temperature-correction procedure is geometry-independent, and so these methods readily generalize to 2- and 3-D problems.

For an atmosphere in differential motion, the energy-conservation constraint yields a radiative flux that is rigorously divergence-free in every local matter frame. Determining the temperature stratification by bringing matter into thermal equilibrium with such a radiation field - i.e., by imposing radiative equilibrium in the co-moving frame - is an excellent approximation if the local cooling time scale is short compared to the local expansion time scale. This condition is well satisfied for the spectrum-forming layers of supernovae (SNe) and of hot star winds (Klein & Castor 1978).

The constraint that the energy packets be indivisible is advantageous from the point of view of coding simplicity. The interaction histories of the packets are then followed one-by-one as they propagate through the computational domain, with there being no necessity to return to any of a packet’s interactions in order to continue or complete that interaction. This is to be contrasted with a Monte Carlo code that directly simulates physical processes by taking its quanta to be a sampling of the individual photons. Absorption of a Monte Carlo quantum is then often followed by the emission of several quanta as an atom cascades back to its ground state. Multiple returns to this interaction are then necessary in order to follow the subsequent paths of each of these cascade quanta. The resulting coding complexity is of course compounded by some of these quanta creating further cascades.

Although coding simplicity argues strongly for indivisible packets, a counter argument is the apparent implied need to approximate the treatment of line formation. Thus, in Monte Carlo codes for studying the dynamics of stellar winds (Abbott & Lucy 1985; Lucy & Abbott 1993) or for synthesizing the spectra of SNe (Lucy 1987; Mazzali & Lucy 1993), the integrity of the packets could readily be maintained since lines were assumed to form by coherent scattering in the matter frame. But significantly, an improved SN code has recently been described (Lucy 1999b) in which branching into the alternative downward transitions is properly taken into account without sacrificing indivisibility. Accordingly, an obvious question now is whether Monte Carlo techniques can be developed that enforce energy-packet indivisibility and yet do not have to adopt any simplifications with regard to line formation. If this can be achieved, then Monte Carlo codes for general NLTE transfer problems become feasible.

## 2 Macro-atoms

As discussed in Sect. 1, it is common in Monte Carlo transfer codes to quantize radiation into monochromatic energy packets. But matter is not quantized, neither naturally into individual atoms nor artificially into parcels of matter. Instead, the continuum description of matter is retained, with macroscopic absorption and scattering coefficients governing the interaction histories of the energy packets.

Nevertheless, it now proves useful to imagine that matter is quantized into macro-atoms whose properties are such that their interactions with energy packets asymptotically reproduce the emissivity of a gas in statistical equilibrium. But these macro-atoms, unlike energy packets, do not explicitly appear in the Monte Carlo code. As conceptual constructs, they facilitate the derivation and implementation of the Monte Carlo transition probabilities that allow in an accurate treatment of line formation.

The general properties of macro-atoms are as follows:

1) Each macro-atom has discrete internal states in one-to-one correspondence with the energy levels of the atomic species being represented.

2) An inactive macro-atom can be activated to one of its internal states $`i`$ by absorbing a packet of kinetic energy or a packet of radiant energy of an appropriate co-moving frequency.

3) An active macro-atom can undergo an internal transition from state $`i`$ to any other state $`j`$ without absorbing or emitting an energy packet.

4) An active macro-atom becomes inactive by emitting a packet of kinetic energy or a packet of radiant energy of an appropriate co-moving frequency.

5) The de-activating packet has the same energy in the macro-atom’s frame as the original activating packet.

Figure 1 illustrates these general rules. An inactive macro-atom, with internal states shown schematically, encounters a packet of energy $`ϵ_0`$ and is activated to one of these states. The active macro-atom then undergoes two internal transitions before de-activating itself by emitting a packet of energy $`ϵ_0`$.

Subsequently, energy packets will in general be referred to as $`e`$-packets but also as $`r`$\- or $`k`$-packets when specifying their contents to be radiant or kinetic energy, respectively.

## 3 Transition probabilities

In Sect. 2, the concept of a macro-atom was introduced by stating some general properties concerning its interaction with $`e`$-packets. The challenge now is to derive explicit rules governing a macro-atom’s activation, its subsequent internal transitions, and its eventual de-activation. Asymptotically, the result of obeying these rules must be the emissivity corresponding to statistical equilibrium.

### 3.1 Energy flow rates

For the moment, we drop the notion of a macro-atom and consider a real atomic species interacting with its environment. Let $`ϵ_i`$ denote the excitation plus ionization energy of level $`i`$ and let $`R_{ij}`$ denote the radiative rate for the transition $`ij`$. The rates per unit volume at which transitions into and out of $`i`$ absorb and emit radiant energy are then

$$\dot{A}_i^R=R_\mathrm{}iϵ_i\mathrm{}and\dot{E}_i^R=R_i\mathrm{}ϵ_i\mathrm{},$$
(1)
respectively, where $`ϵ_i\mathrm{}=h\nu _i\mathrm{}=ϵ_iϵ_{\mathrm{}}`$. Note the summation convention adopted for the suffix $`\mathrm{}`$, which ranges over all levels $`<i`$, including those of lower ions. Similarly, below, the suffix $`u`$ implies summation over all levels $`>i`$, including those of higher ions.

The corresponding rates at which kinetic energy is absorbed from, or contributed to, the thermal pool by transitions to and from level $`i`$ are

$$\dot{A}_i^C=C_\mathrm{}iϵ_i\mathrm{}and\dot{E}_i^C=C_i\mathrm{}ϵ_i\mathrm{},$$
(2)
where $`C_{ij}`$ is the collisional rate per unit volume for the transition $`ij`$.

If we now define the total rate for the transition $`ij`$ to be $`_{ij}=R_{ij}+C_{ij}`$, then the net rate at which level $`i`$ absorbs energy is

$$\dot{A}_i^R+\dot{A}_i^C\dot{E}_i^R\dot{E}_i^C=(_\mathrm{}i_i\mathrm{})(ϵ_iϵ_{\mathrm{}}).$$
(3)
This is an identity that follows directly from the defining Eqs.(1) and (2); it is therefore quite general and does not assume statistical equilibrium.

### 3.2 Statistical equilibrium

We now assume that the level populations $`n_i`$ are in statistical equilibrium. For level $`i`$, this implies that

$$(_\mathrm{}i_i\mathrm{})+(_{ui}_{iu})=0.$$
(4)
A useful alternative representation of statistical equilibrium is obtained by multiplying Eq.(4) by $`ϵ_i`$ and then eliminating the term $`(_\mathrm{}i_i\mathrm{})ϵ_i`$ using Eq.(3). The result can be written in the form

$`\dot{E}_i^R+\dot{E}_i^C+_{iu}ϵ_i+_i\mathrm{}ϵ_{\mathrm{}}`$
$`=\dot{A}_i^R+\dot{A}_i^C+_{ui}ϵ_i+_\mathrm{}iϵ_{\mathrm{}}.`$ (5)
Eq.(4), the conventional equation of statistical equilibrium, balances the rates at which basic atomic processes excite and de-excite level $`i`$. As such, it directly relates to Nature’s quantization of radiation into photons and of matter into atoms. In contrast, Eq.(5), though mathematically equivalent, deals with macroscopic energy flow rates in a finite volume element. These flows can now be quantized into indivisible $`e`$-packets. Moreover, we can think of the volume element as a macro-atom with discrete energy states.

### 3.3 Interpretation

Eq.(5) expresses the fact that in statistical equilibrium the contribution from level $`i`$ to the energy content of unit volume is stationary. In consequence, the net rate at which level $`i`$ gains energy - the right-hand side of Eq.(5) - equals the net rate of loss - the left-hand side.

But the importance here of Eq.(5) lies in the various terms contributing to gains and losses by level $`i`$ and their relevance for constucting transition rules for macro-atoms. The net rate of gain comprises the expected absorption terms $`\dot{A}_i^R`$ and $`\dot{A}_i^C`$ plus the terms $`_{ui}ϵ_i`$ and $`_\mathrm{}iϵ_{\mathrm{}}`$ that clearly represent energy flows into $`i`$ from upper and lower levels. Similarly, the net rate of loss comprises the expected emission terms $`\dot{E}_i^R`$ and $`\dot{E}_i^C`$ plus the terms $`_{iu}ϵ_i`$ and $`_i\mathrm{}ϵ_{\mathrm{}}`$ representing energy flows out of $`i`$ to upper and lower levels.

The above remarks imply definitive values for the energy flows between level $`i`$ and other levels. But this is not true. If Eq.(4) is rewritten as

$$_{iu}+_i\mathrm{}=_{ui}+_\mathrm{}i,$$
(6)
then comparison with Eq.(5) shows immediately that an arbitrary quantity of energy $`ϵ`$ may be added to $`ϵ_i`$ and $`ϵ_{\mathrm{}}`$ without invalidating this equation. But this merely shifts the zero point of the energy scale for excitation and ionization, which we are always free to do. Nevertheless, this freedom implies a corresponding indefiniteness in the energy flow rates between levels.

### 3.4 Stochastic interpretation

Notwithstanding this indefiniteness, we now interpret Eq.(5) in terms of macro-atoms absorbing and emitting $`e`$-packets or undergoing transitions between internal states. In this interpretation, the terms $`\dot{A}_i^R`$ and $`\dot{A}_i^C`$ obviously represent the activation of macro-atoms to state $`i`$ due to the absorption of $`r`$-packets and of $`k`$-packets, respectively.

Now consider an ensemble of active macro-atoms in state $`i`$. For this ensemble to reproduce the behaviour of the real system, the relative numbers of the macro-atoms that subsequently de-activate with the emission an $`r`$\- or $`k`$-packet or which make a transition to another internal state must be proportional to the relative values of the corresponding terms on the left-hand side of Eq.(5). Accordingly, for an individual macro-atom in state $`i`$, the probabilities that it de-activates with the emission of an $`r`$-packet or a $`k`$-packet are

$$p_i^R=\dot{E}_i^R/D_iandp_i^C=\dot{E}_i^C/D_i,$$
(7)
where

$$D_i=\dot{E}_i^R+\dot{E}_i^C+_{iu}ϵ_i+_i\mathrm{}ϵ_{\mathrm{}}=(_i\mathrm{}+_{iu})ϵ_i.$$
(8)
Similarly, the probabilities that it makes an internal transition to particular upper or lower states are

$$p_{iu}=_{iu}ϵ_i/D_iandp_i\mathrm{}=_i\mathrm{}ϵ_{\mathrm{}}/D_i.$$
(9)
Unlike transition probabilities for real atoms, these analogues for macro-atoms depend on ambient conditions. Consequently, in the course of a NLTE calculation, they are iterated on just as are Eddington factors in various other radiative tranfer schemes (Auer & Mihalas 1970; Hummer & Rybicki 1971). Moreover, as with Eddington factors, the Monte Carlo transition probabilities are dimensionless ratios that are likely to converge faster than do their dimensional numerators and denominators.

### 3.5 Excitation equilibrium

When Eq.(5) is summed over all energy levels, the energy flows between different levels cancel, giving

$$\underset{i}{}(\dot{A}_i^R+\dot{A}_i^C)=\underset{i}{}(\dot{E}_i^R+\dot{E}_i^C)$$
(10)
Thus, in statistical equilibrium, the energy stored in the form of excitation and ionization is stationary. For the macro-atoms, this is obeyed rigorously by each activation - de-activation event since the emitted packet’s energy equals that of the absorbed packet - see Figure 1.

## 4 Alternative formulations

Monte Carlo transition probabilities have been defined in Sect. 3, but their non-negativity was not established. Of concern in this regard is stimulated emission when level populations are inverted. However, in anticipation of this issue, radiative rates were introduced without specifying whether stimulated emission contributes positively to $`R_{ij}`$ or negatively to $`R_{ji}`$. We now exploit this flexibility in order to avoid negative probabilities.

### 4.1 General case

In the general case, inverted level populations may occur - i.e., $`g_jn_i>g_in_j`$ for some $`i>j`$ .

#### 4.1.1 Definitions of rates

In order to prevent the probabilities becoming negative when levels invert, stimulated emissions must be added to spontaneous emissions and not treated as negative absorptions. Accordingly, for bound-bound (b-b) transitions, the radiative rates per unit volume are defined to be

$$R_{ij}=(A_{ij}+B_{ij}\overline{J}_{ij}^e)n_iandR_{ji}=B_{ji}\overline{J}_{ji}^an_j,$$
(11)
where $`\overline{J}_{ij}^e`$ and $`\overline{J}_{ji}^a`$ are the mean intensities averaged over the line’s emission and absorption profiles \- see Mihalas (1978, p78). Similarly, for free-bound (f-b) and bound-free (b-f) transitions, we define

$$R_{\kappa i}=(\alpha _i^{sp}+\alpha _i^{st})n_\kappa n_eandR_{i\kappa }=\gamma _in_i.$$
(12)
Here $`\alpha _i^{sp}`$ and $`\alpha _i^{st}`$ are the rate coefficients for spontaneous and stimulated recombinations to level $`i`$, and $`\gamma _i`$ is the uncorrected rate coefficient for photoionizations from level $`i`$. Each of these three quantities can be expressed as an integral over frequency involving the b-f absorption coefficient for an atom excited to level $`i`$ \- see Mihalas (1978, pp130-131).

For collisions, a population inversion gives a negative rate if de-excitations are treated as negative excitations. This is avoided by defining

$$C_{ij}=q_{ij}n_in_eandC_{ji}=q_{ji}n_jn_e.$$
(13)
With these expressions for the radiative and collisional rates, the probabilities defined by Eqs.(7) and (9) are non-negative provided only that the $`ϵ_i`$’s are non-negative. This latter condition is satisfied with the standard convention that the ground state of the neutral atom has zero excitation energy.

#### 4.1.2 Absorption of packets

Because $`R_\mathrm{}i`$ and therefore $`\dot{A}_i^R`$ are here defined without correcting for stimulated emission, the macroscopic line- and continuum-absorption coefficients that determine the flight paths of $`r`$-packets must also be defined without this correction. This ensures a positive absorption coefficient even for a transition with a population inversion.

#### 4.1.3 Emission of packets

If the Monte Carlo transition probabilities result in a macro-atom de-activating radiatively from state $`i`$, the next step is to determine the frequency of the photons comprising the emitted $`r`$-packet. First we suppose that $`i`$ corresponds to a bound level.

Because $`R_i\mathrm{}`$ and therefore $`\dot{E}_i^R`$ here include stimulated emission, the process that radiatively de-activates the macro-atom may be either a spontaneous or a stimulated emission. The ratio of the probabilities of these alternatives is $`q=\dot{E}_i^{sp}/\dot{E}_i^{st}`$, where

$$\dot{E}_i^{sp}=A_i\mathrm{}n_iϵ_i\mathrm{}=\dot{E}_i\mathrm{}^{sp}and\dot{E}_i^{st}=B_i\mathrm{}\overline{J}_i\mathrm{}^en_iϵ_i\mathrm{}=\dot{E}_i\mathrm{}^{st}$$
(14)
are the contributions to $`\dot{E}_i^R`$ from spontaneous and stimulated emissions. Knowing $`q`$, we can choose between the two alternatives with a standard Monte Carlo procedure. Thus, if $`x`$ is a random number from the interval $`(0,1)`$, we select spontaneous emission if $`x<q/(1+q)`$ and stimulated otherwise.

Having thus decided the emission process, we must next choose a downward transition. For spontaneous line emission, the transition $`ij`$ is selected with probability $`\dot{E}_{ij}^{sp}/\dot{E}_i^{sp}`$. For stimulated emission, on the other hand, the selection probability is $`\dot{E}_{ij}^{st}/\dot{E}_i^{st}`$.

With the transition thus determined, the frequency $`\nu `$ of the $`r`$-packet is selected by randomly sampling the line’s emission profile $`\varphi _\nu ^e`$. Thus, if $`x`$ again denotes a random number from $`(0,1)`$, then $`\nu `$ is determined by the equation

$$_0^\nu \varphi _\nu ^e𝑑\nu =x.$$
(15)
This equation can of course always be solved numerically for $`\nu `$. However, elegant and efficient procedures for sampling standard profiles are available (Lee 1974a,b).

Now we consider a macro-atom that de-activates from a continuum state $`\kappa `$. In this case, the probabilities of spontaneous and stimulated emission are in the ratio $`\dot{E}_\kappa ^{sp}:\dot{E}_\kappa ^{st}`$, where

$$\dot{E}_\kappa ^{sp}=\alpha _i\mathrm{}^{sp}n_\kappa ϵ_\kappa \mathrm{}and\dot{E}_\kappa ^{st}=\alpha _i\mathrm{}^{st}n_\kappa ϵ_\kappa \mathrm{}$$
(16)
are the contributions to $`\dot{E}_\kappa ^R`$ from spontaneous and stimulated emissions. Thus $`\nu `$ is selected by first deciding between spontaneous and stimulated emission and then randomly sampling the energy distribution of the chosen process’s recombination continua.

#### 4.1.4 Direction of propagation

If the above selection procedure rules that an $`r`$-packet is emitted spontaneously, then a new direction of propagation is chosen in accordance with this process’s isotropic emission. On the other hand, for stimulated emission, the new direction of propagation is that of the stimulating photon. Thus, the new direction will be in solid angle $`d\omega `$ at $`\theta ,\varphi `$ with probability $`d\omega /4\pi \times I_\nu (\theta ,\varphi )/J_\nu `$, where $`\nu `$ is the frequency of the emitted $`r`$-packet. Accordingly, a Monte Carlo code that treats stimulated emission separately must store a complete description of the radiation field - i.e., $`I_\nu (\theta ,\varphi )`$.

### 4.2 Standard case

For problems where population inversions are not anticipated, we can usefully make the traditional assumption that lines have identical emission and absorption profiles and treat stimulated emissions as negative absorptions - see Mihalas (1978, p78).

#### 4.2.1 Definitions of rates

The radiative rates for b-b transitions are then

$$R_{ij}=A_{ij}n_iandR_{ji}=(B_{ji}n_jB_{ij}n_i)\overline{J}_{ji}.$$
(17)
Similarly, for f-b and b-f transitions, we define

$$R_{\kappa i}=\alpha _i^{sp}n_\kappa n_eandR_{i\kappa }=\gamma _i^{corr}n_i,$$
(18)
where the photionization coefficient is now corrected for stimulated recombinations.

For collisions, the absence of population inversions allows us to treat de-excitations as negative excitations without the risk that Eqs.(7) and (9) will give negative probabilities. Accordingly, we now define

$$C_{ij}=0andC_{ji}=(q_{ji}n_jq_{ij}n_i)n_e.$$
(19)
This then implies that $`\dot{E}_i^C`$ and therefore also $`p_i^C=0`$ for all $`i`$. Energy transfer from the radiation field to the thermal pool then occurs explicitly only via f-f absorptions.

#### 4.2.2 Absorption of packets

Because $`R_\mathrm{}u`$ and therefore $`\dot{A}_i^R`$ are here defined with the correction for stimulated emission included, the macroscopic line- and continuum-absorption coefficients must also include this correction. In the posited absence of population inversions, these absorption coefficients are positive.

#### 4.2.3 Emission of packets

Because $`R_i\mathrm{}`$ and therefore $`\dot{E}_i^R`$ now exclude stimulated emission, the process that radiatively de-activates a macro-atom is always a spontaneous emission. If $`i`$ is a bound state, the frequency $`\nu `$ of the emitted $`r`$-packet is then decided as follows: the transition $`ij`$ is selected with probability $`A_{ij}n_iϵ_{ij}/\dot{E}_i^R`$, and then $`\nu `$ is selected by randomly sampling this transition’s emission profile, as in Sect. 4.1.3.

For de-activation from a continuum state, $`\nu `$ is selected by randomly sampling the energy distribution of the spontaneous recombination continua.

#### 4.2.4 Direction of propagation

Because the de-activating process is in this case spontaneous emission, the new direction of propagation is selected according to isotropic emission. Thus, we now do not need to store $`I_\nu (\theta ,\varphi )`$. In fact, from the Monte Carlo radiation field generated at one iteration, we only require the mean intensities $`J_\nu `$. These allow us to compute transition probabilities from Eqs.(7) and (9) for use during the next iteration.

### 4.3 Large velocity gradients

The procedures described in Sects. 4.1 and 4.2 apply to both static and moving media. But for some important problems involving moving media, a substantial speeding up of the calculation with negligible loss of accuracy is possible by applying Sobolev’s theory of line formation. In doing so, we take advantage of a small dimensionless quantity - the ratio of a line’s Doppler width to the typical flow velocity, which implies an essentially constant velocity gradient over the zone in which a given pencil of radiation interacts with a particular line. The Monte Carlo codes for hot star winds and SNe cited in Sect. 1 treat line formation in the Sobolev approximation.

The simplest case of this kind is that of homologous spherical expansion, as is commonly assumed for SNe. This case will be treated here since it will be used in the test calculations of Sect. 5. But generalization to a spherically-symmetric stellar wind is readily carried out by referring to Castor & Klein (1978). We also assume no population inversions and so treat stimulated emissions as negative absorptions, as in Sect. 4.2.

#### 4.3.1 Definitions of rates

The radiative rates for b-b transitions are then

$$R_{ij}=A_{ij}\beta _{ji}n_iandR_{ji}=(B_{ji}n_jB_{ij}n_i)\beta _{ji}J_{ji}^b.$$
(20)
Here $`J_{ji}^b`$ is the mean intensity at the far blue wing of the transition $`ji`$, and $`\beta _{ji}`$ is the Sobolev escape probability for this transition, given by

$$\beta _{ji}=\frac{1}{\tau _{ji}}[1\mathrm{exp}(\tau _{ji})],$$
(21)
where $`\tau _{ji}`$, the transition’s Sobolev optical depth, is

$$\tau _{ji}=(B_{ji}n_jB_{ij}n_i)\frac{hct_E}{4\pi },$$
(22)
with $`t_E`$ being the elapsed time since the SN exploded. For f-b and b-f transitions, the rates are as in Eq.(18). For collisions, the rates are as in Eq.(19).

#### 4.3.2 Absorption of packets

The absorption of $`r`$-packets by lines is determined by the Sobolev optical depths given by Eq.(22). Absorption of an $`r`$-packet to the continuum is determined by the conventional macroscopic absorption coefficient corrected for stimulated emission.

#### 4.3.3 Emission of packets

The frequency of an emitted $`r`$-packet is decided as follows: for de-activation from a bound state $`i`$, the transition $`ij`$ is selected with probability $`A_{ij}\beta _{ji}n_iϵ_{ij}/\dot{E}_i^R`$, where $`\dot{E}_i^R`$ is evaluated with Eq.(1) using the decay rates from Eq.(20), and the emitted packet is assigned frequency $`\nu _{ij}^{}`$ \- i.e., it is in the far red wing of a line whose emission profile is approximated by a delta function. The packet’s next possible b-b transition is therefore with the next line to the redward of $`\nu _{ij}`$ (Abbott & Lucy 1985).

For de-activation from a continuum state, the new frequency is, as in Sect. 4.2.3, selected by randomly sampling the energy distribution of the spontaneous recombination continua.

#### 4.3.4 Direction of propagation

If an $`r`$-packet is emitted from a continuum state, the new direction of propagation is selected according to isotropic emission since the emission in this case is spontaneous. For de-activation from a bound state, the emission is also isotropic since, for homologous expansion, there is no kinematically-preferred direction. This is not true for a stellar wind.

## 5 Convergence tests

The Monte Carlo transition probabilities derived in Sect. 3 are designed to reproduce asymptotically the emissivity of an atomic species whose level populations are in statistical equilibrium. To test this, we now consider one-point problems with specified and fixed ambient conditions. Such tests sensibly precede application to a general NLTE problem, for then the local ambient conditions are everywhere being adjusted iteratively as the global solution is sought.

### 5.1 Fe II

In the initial tests, the Monte Carlo transition probabilities are applied to the model Fe II ion with $`N=394`$ levels used previously (Lucy 1999b) to investigate the accuracy of approximate treatments of line formation in SNe envelopes. The energy levels of the Fe II ion and the f-values for permitted transitions were extracted from the Kurucz–Bell (1995) compilation by M.Lennon (Munich). Einstein A-values for forbidden transitions are from Quinet et al.(1996) and Nussbaumer & Storey (1988). Collision strengths, needed for Sect.5.1.5, are from Zhang & Pradhan (1995) and van Regemorter (1962).

#### 5.1.1 Radiative excitation

In the first Fe II test, we neglect collisional excitations and, as previously (Lucy 1999b), take the ambient radiation field determining the quantities $`J_{ji}^b`$ in Eq.(20) to be $`WB_\nu (T_b)`$ with $`T_b=12500K`$ and dilution factor $`W=0.5`$, corresponding to $`r=R`$. The density parameter is $`n(FeII)=6.6\times 10^7cm^3`$, and the time since explosion is $`t_E=13days`$. With parameters specified, this one-point statistical equilibrium problem - Eq.(4) for $`N1`$ levels plus a normalization constraint - is non-linear in the unknowns $`n_i`$ because the rate coefficients in Eq.(20) depend on the $`n_i`$ through the Sobolev escape probabilities. Fortunately, repeated back substitutions give a highly accurate solution $`n_i^{(x)}`$ in $``$ 10 iterations.

#### 5.1.2 Monte Carlo experiment

With $`n_i^{(x)}`$ determined, the Fe II level emissivites $`\dot{E}_i^R`$ and absorption rates $`\dot{A}_i^R`$ can be computed from Eq.(1). We now test the Monte Carlo transition probabilities by seeing how accurately they reproduce these values $`\dot{E}_i^R`$. Note that it is sufficient to test level emissivities since if these are exact so also are the line emissivities computed as described in Sect. 4.3.3.

In the following Monte Carlo experiment, $`𝒩`$ packets of radiant energy are absorbed and subsequently emitted by a macro-atom representing a macroscopic volume element of Fe II ions in the ambient conditions specified above. The energies of these packets are taken to be equal and given by $`ϵ_0=\dot{A}^R/𝒩`$, where $`\dot{A}^R=_i\dot{A}_i^R`$. The calculation proceeds step-by-step as follows:

1) $`𝒩_i=𝒩\dot{A}_i^R/\dot{A}^R`$ of the packets activate the macro-atom to internal state $`i`$.

2) The transition probabilities $`p_i^R`$, $`p_{iu}`$ and $`p_i\mathrm{}`$ for a macro-atom in state $`i`$ are computed from Eqs.(7) and (9).

3) The transition probabilities sum to one, so each corresponds to a segment $`(x_k,x_{k+1})`$ of the interval (0,1). A particular transition is therefore selected by computing a random number $`x`$ in (0,1) and finding in which segment it falls.

4) If the selected transition is the de-activation of the macro-atom, we update $`\dot{E}_i^{MC}`$ to $`\dot{E}_i^{MC}+ϵ_0`$ and then return to step 3) to process the next activation of state $`i`$, or to step 2) to process the first of the packets that activate the macro-atom to state $`i+1`$.

5) If the selected transition is an internal transition to state $`j`$, then we return to step 2) with $`j`$ replacing $`i`$.

6) When all $`𝒩`$ packets have been processed, the final elements of the vector $`\dot{E}_i^{MC}`$ are the estimates of the level emissivities $`\dot{E}_i^R`$

#### 5.1.3 Results of experiment

As a single measure of the accuracy of the estimated level emissivities, we compute the quantity

$$\delta =\underset{i}{}|\dot{E}_i^{MC}\dot{E}_i^R|/\underset{i}{}\dot{E}_i^R.$$
(23)
This is the mean of the absolute fractional errors of the $`\dot{E}_i^{MC}`$ when weighted by $`\dot{E}_i^R`$.

Figure 2 shows the values of $`\delta `$, expressed as percentage errors, found in a series of trials with $`𝒩`$ increasing from $`10^4`$ to $`10^7`$. The values of $`\delta `$ decrease monotonically with increasing $`𝒩`$, falling to $`0.36`$ percent for $`𝒩=10^7`$. More importantly, the errors accurately follow an $`𝒩^{1/2}`$ line, as expected if the only source of error are the sampling error at step 3) of the Monte Carlo experiments. Accordingly, to the accuracy of these experiments, macro-atoms obeying the transition probabilities derived in Sect. 3 do indeed reproduce the emissivity of a gas in statistical equilibrium.

Also included in Fig.2 are values of $`\delta `$ obtained when the transition probabilities are computed with excitation energies $`ϵ_i`$ increased by 5eV. This is to investigate the consequences of the dependence of the energy flow terms in Eq.(5) - and therefore also of the transition probabilities - on the zero point of the scale of excitation energy. These results also track an $`𝒩^{1/2}`$ line and so indicate that the predicted emissivities are asymptotically independent of the zero point. But since the open circles are marginally higher, there is an indication that increasing the zero point gives slighty less accurate emissivities at a given $`𝒩`$.

In the Monte Carlo codes for hot star winds and SNe cited in Sect. 1, line formation is treated approximately, with either resonant scattering or downward branching being assumed. For both assumptions, $`\dot{E}_i^R=\dot{A}_i^R`$, corresponding to a macro-atom for which de-activation always immediately follows activation - i.e., $`p_i^R=1`$ for all $`i`$. In this case, as indicated on Fig.2, $`\delta =7.28`$ percent. Thus, when the points in Fig.2 drop below this value, the success must be due to the internal, radiationless transitions governed by the probabilities $`p_{iu}`$ and $`p_i\mathrm{}`$.

#### 5.1.4 Distribution of jumps

The above experiments show that despite the formidable complexity of its level structure the Fe II ion’s reprocessing of radiation is accurately simulated by the Monte Carlo transition probabilities. Nevertheless, from a computational standpoint, a remaining concern is how many internal transitions - or jumps - does this require? To answer this, the number of jumps before de-activation was recorded for each absorbed packet in the $`𝒩=10^7`$ trial and used to derive $`N(j)`$, the number of packets requiring $`j`$ jumps.

From $`N(j)`$, we find that the expected number of jumps is $`<j>=2.19`$ and that the probability of immediate de-activation - i.e., zero jumps - is $`P_0=0.425`$. Evidently, fears of numerous, time- consuming internal transitions are ill-founded.

Figure 3 is a logarithmic plot of $`N(j)`$. This reveals a power-law decline with increasing $`j`$ but with alternating deviations indicating that an even number of jumps before de-activation is favoured. A simple model suggests the origin of this curious behaviour. Consider a 3-level atom with $`ϵ_3>ϵ_2>ϵ_1=0`$ and suppose that level 2 is metastable with $`A_{21}=0`$. Because $`B_{12}=0`$, the macro-atom can only be activated to state 3; and because $`A_{21}=B_{21}=0`$, the macro-atom cannot de-activate from state 2. Moreover, since $`ϵ_1=0`$, Eq.(9) gives $`p_{31}=p_{21}=0`$, and so state 1 of the macro-atom cannot be reached. Accordingly, following activation at state 3, the macro-atom de-activates with probability $`p`$ or jumps to state 2 with probability $`1p`$, from whence it returns to state 3 with probability $`p_{23}=1`$. It is now simple to prove that the probabilty of $`j`$ jumps before de-activation is $`P_j=p(1p)^{j/2}`$ if $`j`$ is even, and $`P_j=0`$ if $`j`$ is odd. The Fe II ion’s numerous low-lying metastable levels are presumably playing the role of level 2 and thereby favouring an even number of jumps.

Histograms $`N(j)`$ have also been computed for two other cases. First, the above trial was repeated with the $`ϵ_i`$’s increased by 5eV as in Sect. 5.1.3. This change increases $`<j>`$ \- to 4.54 - as expected since the probabilities $`p_{iu}`$ and $`p_i\mathrm{}`$ are thereby increased and $`p_i^R`$ correspondingly decreased. Evidently, the standard choice of energy-level zero point leads to the most computationally-efficient set of transition probabilities.

In the second case, $`W`$ is decreased from 0.5 to 0.067, corresponding to $`r=2R`$. This change decreases $`<j>`$ \- from 2.19 to 1.29 - as expected given the weakening of the radiative excitation rates.

#### 5.1.5 Collisional excitation

In the above experiment, the emission derives entirely from radiative excitation since collisions were neglected. Now we consider the opposite extreme by setting the ambient radiation field to zero but including collisions.

The only parameters of this test are the electron temperature and density, and these are assigned the values $`T_e=2\times 10^4K`$ and $`N_e=10^8cm^3`$. The resulting statistical equilibrium problem is linear and so solved without iteration. For this solution, accurate values of the level emissivities $`\dot{E}_i^R`$ are again computed from Eq.(1).

The next step is to derive estimates of the level emissivities by repeating the Monte Carlo experiment of Sect.5.1.2. The only changes needed are the following: first, since the solution has population inversions the general formulation of Sect. 4.1 must be adopted to avoid negative probabilities.

Secondly, since a macro-atom is now always activated by a $`k`$-packet, their energies are taken to be $`ϵ_0=\dot{A}^C/𝒩`$, where $`\dot{A}^C=_i\dot{A}_i^C`$. Correspondingly, at step 1) of the experiment, $`𝒩_i=𝒩\dot{A}_i^C/\dot{A}^C`$.

Thirdly, since a macro-atom can now de-activate by emitting either an $`r`$\- or a $`k`$-packet, only the former results in an updating of $`\dot{E}_i^{MC}`$. The emission of a $`k`$-packet represents the return of energy $`ϵ_0`$ to the therrmal pool.

Apart from these changes, the convergence experiment proceeds as in Sects. 5.1.2 and 5.1.3. The result is a plot similar to Fig.2, but with $`\delta =0.19`$ percent for $`𝒩=10^7`$. Evidently, the Monte Carlo transition probabilities are equally applicable to problems where collisional excitation is a source of emission.

### 5.2 Hydrogen

Although the Fe II experiments demonstrate the validity of the Monte Carlo transition probabilities, a test including b-f and f-b transitions is of interest. Accordingly, a convergence experiment at one point in a SN’s envelope has also been carried out for a 15-level model of the H atom, with level 15 being the continuum $`\kappa `$. The 14 bound levels correspond to principal quantum numbers $`n=114`$, with each level having consolidated statistical weight $`g=2n^2`$.

As for Fe II, the ambient radiation field incident on the blue wings of the b-b transitions is $`WB_\nu (T_b)`$, but now with $`T_b=6000K`$ and $`W=0.067`$. However, beyond the Lyman limit, we assume zero intensity, so that photoionizations occur only from excited states. Correspondingly, recombinations to $`n=1`$ are excluded on the assumption of immediate photoionization. Collisional excitations and ionizations are neglected. The density parameter is $`N(H)=1.88\times 10^9cm^3`$, the electron temperature $`T_e=4800K`$, and the time since explosion $`t_E=10days`$. With parameters specified, this non-linear statistical equilibrium problem can also be solved with repeated back substitutions, giving a highly accurate solution $`n_i^{(x)}`$ in $`30`$ iterations.

With $`n_i^{(x)}`$ determined, Monte Carlo experiments as described in Sect. 5.1.2 were carried out to test if level emissivities are also recovered in this case. In Fig.4, two such trials, with $`𝒩=10^4`$ and $`10^5`$, are compared with the exact solution. The results show that excellent agreement is achieved for $`𝒩=10^5`$. Note in particular the success with $`\dot{E}_\kappa ^R`$, which is the rate of ionization energy loss due to recombinations, and with $`\dot{E}_2^R`$, whose very low value is due to the strong trapping of $`L\alpha `$ photons.

### 5.3 Alternative test of convergence

Thus far, a Monte Carlo procedure has been used to validate the transition probabilities developed in Sect.3. This has the advantage of following closely and therefore illustrating their use in realistic NLTE calculations. But for feasible values of $`𝒩`$, sampling errors limit the accuracy of such tests.

In order to test to higher precision, approximate level emissivities $`\dot{E}_i^{(m)}`$ can be computed recursively according to the following scheme:

$$\dot{E}_i^{(m)}=p_i^R\underset{r=1}{\overset{m}{}}G_i^{(r)},$$
(24)
where $`p_i^R`$ is the radiative de-activation probability from Eq.(7) and $`G_i^{(r)}`$ is the increment at cycle $`r`$ to the summation approximating the rate at which level $`i`$ gains energy - i.e. the right-hand side of Eq.(5). This increment is derived from the previous increment by applying the transition probabilities from Eq.(9). Thus

$$G_i^{(r)}=\underset{j}{}p_{ji}G_j^{(r1)},$$
(25)
and the recursion cycles are initiated by setting

$$G_i^{(1)}=\dot{A}_i^R+\dot{A}_i^C.$$
(26)
This procedure is now applied to the Fe II test problem of Sect. 5.1.1. As with that experiment, the accuracy of the vectors $`\dot{E}_i^{(m)}`$ are measured by computing $`\delta `$ defined by Eq.(23). For $`m=17`$, $`\delta `$ drops below the value 0.36 percent found in Sect.5.1.3 with $`𝒩=10^7`$ \- see Fig.2. As the recursion procedure continues further, $`\delta `$ decreases monotonically until at $`m60`$ it drops to a value of $`10^8`$, at which point machine precision or accumulated roundoff errors halt further progress. This test clearly confirms and strengthens the earlier tests of the Monte Carlo transition probabilities.

## 6 Sensitivity tests

The experiments of Sect.5 demonstrate that, when computed with the exact level populations $`n_i^{(x)}`$, the Monte Carlo transition probabilities applied to indivisible $`e`$-packets reproduce the exact level emissivities as $`𝒩\mathrm{}`$. But this success, though necessary, does not of itself imply that the technique will be successful when applied to NLTE problems. For example, if the Monte Carlo emissivities were to undergo large changes in response to small changes in $`n_i`$, then we would reasonably suspect that the iterations inevitably required for a NLTE problem would converge very slowly - or might even diverge. On the other hand, if the emissivities are insensitive to changes in $`n_i`$, then the prospects for successful applications are excellent.

### 6.1 Fe II emissivities

This crucial question of sensitivity can be investigated by repeating the calculations of Fe II emissivities reported in Sect.5.1, but with $`n_i`$ perturbed away from $`n_i^{(x)}`$. A convenient way of doing this is to replace $`n_i^{(x)}`$ by the Boltzmann distribution at excitation temperature $`T_{ex}`$. Then, for given $`T_{ex}`$, the corresponding level emissivities $`\dot{E}_i^{MC}`$ are obtained from a Monte Carlo trial with $`𝒩=5\times 10^6`$ packets, and so are negligibly affected by sampling errors (cf. Fig.2) .

Now, for the given $`T_{ex}`$, we can also compute $`\dot{E}_i^R`$, the level emissivities predicted by the fundamental formulae - Eqs.(1) and (20) in this case. This represents the standard approach to NLTE transfer problems whereby the radiation field is computed from the Radiative Transfer Eq. (RTE) with emissivity coefficients evaluated using the current estimates of $`n_i`$. Thus by comparing these two emissivity estimates $`\dot{E}_i^{MC}`$ and $`\dot{E}_i^R`$, we can see whether this Monte Carlo technique is potentally capable of yielding a superior estimate of the radiation field.

In Fig.5, the quantities $`\dot{E}_i^{MC}`$ and $`\dot{E}_i^R`$ obtained for $`T_{ex}=12500K`$ are plotted against $`\dot{E}_i^{(x)}`$, the exact statistical equilibrium level emissivities \- i.e., the values corresponding to $`n_i^{(x)}`$.

Remarkably, Figure 5 shows that the Monte Carlo emissivities are far less sensitive to the departure of $`n_i`$ from $`n_i^{(x)}`$ than are the emissivities computed directly from the fundamental formula. For the most part, the $`\dot{E}_i^{MC}`$ are in error by $`<0.1`$ dex, with little evidence of bias, while the $`\dot{E}_i^R`$ are systematically offset by $`+0.3`$ dex.

To investigate whether this insensitivity is characteristic of the Monte Carlo procedure, the above test is now repeated with $`T_{ex}`$ ranging from $`7500K`$ to $`20000K`$ and the resulting mean errors defined by Eq.(23) plotted in Fig.6. We see that $`\dot{E}_i^R`$ gives reasonably accurate emissivities only in the immediate neighbourhood of the minimum at $`T_{ex}11250K`$. On the other hand, the values $`\dot{E}_i^{MC}`$ are accurate to $`0.1`$ dex across the entire range.

The causes of these astonishing differences in sensitivity are of considerable interest. For $`\dot{E}_i^R`$, the strong sensitivity to $`T_{ex}`$ is readily understood. Because the sum $`R_i\mathrm{}ϵ_i\mathrm{}n_i`$, an error in the population of the emitting level translates directly into an error in $`\dot{E}_i^R`$.

Now consider $`\dot{E}_i^{MC}`$. This quantity is determined by the rate at which active macro-atoms reach state $`i`$, and this happens by direct absorptions of packets into this state or by transitions from other states. Either way, the accuracy of the source vectors $`\dot{A}_i^R`$ and $`\dot{A}_i^C`$ is clearly fundamental to the accuracy of the vector $`\dot{E}_i^{MC}`$. But the dominant contributors to the elements of these source vectors - see Eqs.(1) and (2) - are transitions from the ground state and from low-lying metastable levels, and the estimated populations of these levels are unlikely to be seriously in error. In particular, with an assumed Boltzmann distribution over excited states, the $`n_i`$ of these low levels is insensitive to $`T_{ex}`$ and do not differ much from $`n_i^{(x)}`$. In contrast, the populations of high levels are quite likely to be badly estimated and are acutely sensitive to $`T_{ex}`$.

### 6.2 Comments

Another way of appreciating the differences in these approaches to calculating emissivities is as follows. The Monte Carlo procedure applies only to a state of statistical equilibrium and, as such, constrains every level’s emissivity to be consistent with the rates of processes populating that level. In contrast, the fundamental emissivity formula applies also to states out of statistical equilibrium and so takes no account of whether the levels’ populations can be maintained. Accordingly, with this Monte Carlo technique, the principle of statistical equilibrium is incorporated (approximately) as the radiation field is being calculated. On the other hand, when emissivities are computed from the fundamental formula, any consideration of statistical equilibrium is effectively being deferred until the updated radiation field has been determined.

The likely beneficial impact of this insensitivity on the iterations needed to derive NLTE solutions is worth stressing. With the conventional RTE approach, an erroneously overpopulated upper level $`i`$ pollutes the radiation field with spurious line photons at frequencies $`\nu _{ij}(j<i)`$, and these are sources of excitation for level $`i`$ when level populations are next solved for. Similarly, an erroneously overpopulated upper ion pollutes the radiation field with recombination photons that are subsequent sources of photoionization for the lower ion. To some degree, therefore, such errors are self-perpetuating and so are not rapidly eliminated. This persistency contributes to the slow convergence typical of NLTE codes. In contrast, with the Monte Carlo approach, this pollution does not happen and so - for sufficiently large $`𝒩`$ \- a high quality radiation field is obtained immediately provided that the initial populations of the low-lying levels are estimated sensibly.

## 7 Implementation

The Monte Carlo transition probabilities allow statistical equilibrium to be incorporated into the calculation of radiation fields for NLTE problems. Moreover, this is achieved without imposing the constraint of radiative equilibrium. Accordingly, in principle at least, the technique applies equally to problems with non-radiative heating, such as stellar chromospheres.

### 7.1 Radiative equilibrium

In the absence of non-radiative heating, a NLTE transfer problem must be solved subject to the constraint of radiative equilibrium. The incorporation of this additional constraint into the macro-atom formalism is readily understood. First suppose that collisional processes are neglected. The absorbed and the emitted $`e`$-packets are then always $`r`$-packets and they have identical energies \- see Fig.1. Thus, the constraint of radiative equilibrium is obeyed rigorously since it holds exactly for every activation - de-activation event, all of which are of the form $`r𝒜^{}r`$, where $`𝒜^{}`$ denotes an active macro-atom. Note also that since active macro-atoms do not appear spontaneously within the computational domain (D), every Monte Carlo quantum’s interaction history starts and ends as an $`r`$-packet crossing a boundary of D.

Now suppose that collisions are included. In this case, a macro-atom activated by an $`r`$-packet can de-activate itself by emitting a $`k`$-packet, so that radiative equilibrium no longer holds exactly for each individual activation - de-activation event. However, the emitted $`k`$-packet is re-absorbed in situ by another macro-atom and thereby (eventually) converted into an $`r`$-packet. Since this has the same energy as the original $`r`$-packet, radiative equilibrium holds for every sequence of in situ events that starts with the absorption of an $`r`$-packet and ends with the next emission of an $`r`$-packet. A typical in situ sequence is $`r𝒜^{}k𝒜^{}k𝒜^{}r`$. If such sequences are abbreviated as $`r[𝒜^{}]r`$, we see that the inclusion of collisions has not fundamentally changed the procedure and that radiative equilibrium is still rigorously obeyed.

### 7.2 Non-radiative heating

In the presence of non-radiative heating, the NLTE problem is not subject to the additional constraint of radiative equilibrium. Statistical equilibrium is incorporated with the macro-atom formalism as before, and the challenge now is to incorporate the creation of radiant energy within D due to the additional heating. This is accomplished by allowing for the spontaneous and random appearance within D of active macro-atoms with their number, locations and internal states $`i`$ all controlled by the collision source vector $`\dot{A}_i^C`$ \- cf. Sect.5.1.5. Note that because this sampling of $`\dot{A}_i^C`$ takes full account of the collisional creation of excitation, the emission of a $`k`$-packet is not now followed by its in situ re-absorption; instead, the interaction history of that Monte Carlo quantum then ends and its energy is added to the thermal pool (cf. Sect.5.1.5.). The radiation field generated by this procedure is not divergence-free but reflects the collisional creation of radiant energy due to an elevated temperature profile maintained by the non-radiative heating.

## 8 Conclusion

The limited aim of this paper has been to see if Monte Carlo transfer codes whose quanta are indestructable energy packets can be constructed without resorting to simplified treatments of line formation. To this end, the concept of a macro-atom has been introduced and rules established governing its activation and de-activation as well as its transitions between internal states. These rules - the Monte Carlo transition probabilities - have been derived by demanding that the macro-atom’s emission of $`r`$-packets asymptotically reproduces the local emissivity of a gas in statistical equilibrium; and these rules’ validity has been confirmed with one-point test problems.

Evidently, the next step is to implement these transition probabilities in a code to solve a realistic NLTE problem for a stratified medium and thus to investigate the practicality of this technique for problems of current interest. In a companion paper, a Monte Carlo NLTE code treating the formation of H lines in a Type II SN envelope will be described and used to illustrate the convergence behaviour of iterations to obtain both the level populations and the temperature stratification.

###### Acknowledgements.

I am grateful to C.Jordan for her interest in the potential application of this technique to stellar chromospheres as well as for a constructive referee’s report.


text:
# On Independent Representation

(Date: October 2005)

## Abstract.

The paper has been withdrawn by the author.