text:
# XMM-Newton observations of OY Car III: OM light curve modelling, X-ray timing and spectral studies

## 1 Introduction

Dwarf novae are interacting stellar binary systems in which material gets transfered from a red dwarf secondary star onto a white dwarf via Roche lobe overflow. In the absence of a strong magnetic field this material forms an accretion disc around the white dwarf. Dwarf novae show outbursts which occur on week to month timescales. For systems which are viewed close to the plane of the binary system, eclipses of the white dwarf by the secondary star are seen every orbital period. These systems are important since they allow a detailed study of the accretion disc as the secondary gradually covers then uncovers the disc.

One such eclipsing dwarf nova is OY Car. It has a binary orbital period of 1.51 hrs and has been extensively studied at various wavelengths. Recently, it has been observed at X-ray and optical/UV energies using XMM-Newton. Ramsay et al (2001a) found that there was a quasi-stable modulation of the X-rays at $``$2240 sec, which was most prominent at lower energies. They speculated that this was the due to disc material orbiting above the binary orbital plane or due to the spin period of the white dwarf. In the latter scenario, this would represent the first time that the spin period of the white dwarf in a dwarf nova had been detected.

The eclipse profile shows a sharp drop in X-rays and in the UV band. In X-rays the flux remains at a very low, but significant level after this rapid (20–30 sec) drop. In the UV and B bands, the flux decreases at a more gradual rate after the rapid drop since the accretion disc is still visible after the eclipse of the white dwarf. Wheatley & West (2003) further analysed the XMM-Newton data and found that the X-ray eclipse lasted 30$`\pm `$3sec and concluded the X-ray emission originated at the poles of the white dwarf.

In contrast, in this paper we investigate the nature of the 2240 sec modulation in more detail by revisiting the XMM-Newton data. In particular we have carried out phased resolved spectroscopy over this modulation period. Further, we apply a new eclipse mapping technique to model the spatial distribution of the accretion disc using the OM XMM-Newton data.

## 2 Observations

XMM-Newton was launched in Dec 1999 and has the largest effective area of any imaging X-ray telescope. It has 3 medium spectral resolution CCD type cameras on-board: two EPIC MOS detectors and one EPIC pn detector. Since the EPIC pn detector has the higher sensitivity we only discuss data taken using this detector (Strüder et al 2001). In addition, it has an 0.3m optical/UV imaging telescope which enables simultaneous X-ray optical/UV observations (Mason et al 2001).

OY Car was observed during the performance verification phase of the mission on June 29-30th 2000. The details of the observation are reported in Ramsay et al (2001a,b). However, the main points to note are that these observations were made some 4 days after an optical outburst. The system had already faded by $`V`$4 mag from its peak, but it was still relatively bright in X-rays ($``$1.1 ct/s in the EPIC pn). The exposure in the EPIC pn detector was 48ksec, it was configured in full frame mode and the medium filter was used. The optical data were obtained in the UVW1 filter (effective wavelength 2910 Å). The data were processed using the XMM-Newton Science Analysis System v5.3.3. Additional optical data in B band was obtained in August 2000 when it was at a similar brightness to the June observation (Ramsay et al 2001a).

## 3 X-ray Light Curve Analysis

Ramsay et al (2001a) found evidence for a quasi-stable X-ray modulation on a period of $``$2240 sec, which was strongest in soft X-rays. This modulation was observed for around 1/3 of the duration of the observation and it was speculated that this period was related to the spin period of the white dwarf.

We have re-examined the light curves of OY Car: in particular we perform a timing analysis of that portion of the light curve in which the 2240 sec modulation is strongest (HJD 2451725.6548–2451725.8622) in 4 energy bands (0.15–10.0, 0.15–1.0, 1–4 and 4–10keV). (We removed the time of the eclipse from the light curves). We show the amplitude spectrum of the 0.15–1.0keV light curve in Figure 1 (top panel). The most prominent peak has a period of 2193 sec: within the errors, this period is consistent with the 2240 sec period found by Ramsay et al (2001a). However, there is a second prominent peak at 3510 sec: this peak is still present if we pre-whiten the data on the 2193 sec period. We performed the same analysis in the 2 more restricted energy bands (Figure 2): we find that both the 2193 sec and 3510 sec periods are more prominent at softer energies. A peak was seen at 3510 sec in the amplitude spectrum shown in Ramsay et al (2001a), however, it is more prominent in this new study since we have restricted the analysis to the time interval in which the quasi-stable modulation is seen. In order to demonstrate the prominence of both of these periods, we show the data folded on either of the two periods, after they have been prewhitened using the other period (Figure 3).

To determine the significance of these two periods we randomly populated the time bins with the observed count rates. We then generated 10000 power density spectra: the 95%, 99% and 99.9% confidence limits from these are shown together with the 0.1-1.0 keV power spectrum in Figure 1(top panel). It is clear that both the 2193 and 3510 sec periods are highly significant. To determine the errors on these periods we used the Cash statistic (Cash 1979), which implies: 2193$`{}_{35}{}^{}{}_{}{}^{+17}`$ and 3510$`{}_{104}{}^{}{}_{}{}^{+25}`$ sec. There are three other periods which are significant at $`>99.9`$% level: 1248$`{}_{15}{}^{}{}_{}{}^{+17}`$ sec, 1574$`{}_{21}{}^{}{}_{}{}^{+18}`$ sec and 4844$`{}_{104}{}^{}{}_{}{}^{+120}`$ sec.

It is possible that some of these less significant periods could be produced as a result of aliasing between the most prominent periods (2193 sec and 3510 sec) and the window function or as a result of these two periods having a non-sinusoidal pulse shape. In order to investigate this further we have simulated a dataset with the same time points as the original data. However, this fake data only contains the two strongest signals (2193 sec and 3510 sec). The shape of these signals is obtained from the real data by using a simultaneous two-period fourier fit, where each of the two pulse shapes are represented with a second order fourier series. As a result we find that at least the 1574s period could be produced by the effects mentioned above, Figure 1 (lower panel).

We discuss the possible origin of all of these periods in §6.1.

## 4 X-ray Spectral Analysis

Taking the same time interval that the quasi-stable behaviour is present as in the previous section, we phased events on the 2193 and 3510 sec quasi-stable periods reported in the previous section. Initially we extracted 6 spectra covering different 2193 sec phase intervals. In extracting these spectra we used apertures $`30^{^{}}`$ in radius centered on OY Car, chosen so that the aperture did not cover more than one CCD. This encompasses $``$90 percent of the integrated PSF (Aschenbach et al 2000). A background spectrum was extracted from the same CCD on which the source was detected, scaled and subtracted from the source spectra.

Since the response of the detectors is not well calibrated at present below $``$0.2keV, energies below this were ignored in the following analysis. The response file epn$`\mathrm{\_}`$ff20$`\mathrm{\_}`$sdY9$`\mathrm{\_}`$medium.rmf was used: this response includes single and double pixel events.

Ramsay et al (2001b) found that they could not well fit the integrated EPIC spectra of OY Car with a single temperature thin plasma model. A multi-temperature model gave a much better fit to the data. In their analysis Ramsay et al (2001b) concluded that a best fit for the integrated spectrum was obtained with a three temperature plasma model with a partially covering absorber. The multi-temperature nature was supported from line H/He-like intensity ratios of Iron and Sulphur and also by the various line species seen in the high resolution RGS spectrum.

We fitted each of the 6 phase resolved spectra separately using the X-ray fitting package XSPEC. Each spectrum was binned so that a minimum of 40 counts were present in each spectral bin. An initial analysis showed that a two temperature fit gave much better fits to the data. In these phase resolved spectra adding a third temperature component did not significantly improve the fit. Further, we find that the temperatures for each of thermal components were consistent in each of the 6 phase resolved spectra (kT1=1.4keV and kT2=10keV). We therefore fixed the temperatures of the thermal plasma models at these temperatures. In each fit the free parameters were the absorption (we used a simple neutral absorption model) and the normalisation. The metal abundance was fixed at solar.

The resulting fits are shown in tabulated form in Table 1 and in graphical form (along with the light curve folded on the 2193 sec period) in Figure 4. The fits are reasonably good in the phase intervals 0.32–0.80, but in the remainder it is less so ($`\chi _\nu ^2`$$``$1.5). In the case of the phase interval 0.0–0.17 the fit was significantly improved when we fix the metal abundance at 1.5 solar ($`\chi _\nu ^2`$=1.35). It is clear from Figure 4 that there is no significant variability in the normalisations of the thermal plasma components. However, there is a variation in the column of the neutral Hydrogen density: it is anti-correlated with the folded intensity light in that it shows a maximum at intensity minimum. We show in Figure 5 the confidence contour plot in the ($`N_H`$, $`T_{low,norm}`$) plane for the spectra covering the phase maximum and minimum: there is a clear separation in the contours.

We show in Figure 6 the X-ray spectra near 6–7keV for spectra taken from the phase maximum and minimum: the profile of the Fe K$`\alpha `$ line does not differ significantly. This further strengthens the argument that the spectral changes over the 2193 sec period are most likely caused by the change in the amount of absorption.

We also folded the data on the 3510 sec quasi-stable modulation: we find a similar result to that found for folding on the 2193 sec. Namely, the temperature and normalisation of the two thermal plasma components do not vary over phase, but the absorption component varies significantly. For brevity we do not show the results in Table or Figure form.

## 5 Modelling the B band and UV data

It is clear from Ramsay et al. (2001a) that the B and UV light curves exhibit clear modulation over the orbital period as well as asymmetric eclipse profiles. In this section we will describe and discuss the modelling of those data in detail.

Our basic assumption is that all of the B and UV modulation is caused by the changing visibility of different parts of (at least) partially optically and physically thick accretion disc and the eclipse of such disc by the secondary star. In our modelling we have taken the idea of ‘fireflies” presented in Hakala, Cropper & Ramsay (2002) and developed it further to suit the modelling of optically thick accretion discs.

Hakala et al. (2002) use a swarm of emitting point sources (fireflies) to model the shape of an accretion stream in eclipsing AM Herculis systems. In their approach the swarm is regularized to favour realistic ‘banana’ shapes for the accretion stream (without fixing the exact location or shape of the stream). They then used a variant of a genetic algorithm to find the optimal accretion stream shape and location.

If we want to apply a similar approach to locating the emitting regions in a disc accreting system, like OY Car, we need to make several modifications to the original code. In fact, the only assumptions that remain the same are the idea of ‘fireflies’ as a numerical model for the emission and the optimisation algorithm. We next discuss our model computation in detail.

### 5.1 The model

We first restrict the possible emission to originate within an accretion disc like volume around the primary. The boundary conditions for this volume are: maximum radius 0.45a, maximum height $`\pm 0.02a`$. In addition we define a somewhat thicker disc like volume of maximum height $`\pm 0.025a`$. This volume is assumed to have a homogenous optical thickness, that remains constant throughout the modelling. The reason for this is to provide estimates for the absorption from the different parts of the disc and thus guarantee an optically thick solution.

To compute the optical light curve for any given distribution of fireflies inside this boundary, the algorithm for the light curve generation procedes as follows. We first loop over all the orbital phases to be fitted and for each phase:

* Select the fireflies that are not eclipsed by the secondary;
* Determine the absorption coefficient for each of the selected fireflies. This is done by computing the length of a path that light from any particular fly has to travel inside a predefined disc volume;
* Sum up the emission from all visible flies at this phase (taking into account the variable absorption for different flies).

The main difference between the ‘firefly approach’ and the ‘normal approach’ is that the firefly approach is (in analogue with hydrodynamics) a particle based solution whilst the ‘normal’ approach relies on fixed grid on which the solution is determined. The reasons behind our choice of the firefly approach here are the following ones. Firstly, one of the main axioms of the standard eclipse mapping routine (Horne 1985) is that the eclipse profile of the disc is solely determined by the changing obscuration of the disc by the secondary. Thus outside the eclipse phases the light curve is constant. This implies that the same algorithm could not be used to fit the rest of the light curve. Normally the out-of-eclipse variation, often seen in quiescent state CV light curves, has to be removed before the fitting, which can cause undesirable bias in the results. Rutten (1998) developed the eclipse mapping by including a 3D shape for the disc. His model could be used for out-of-eclipse light curve modelling, but it has it’s limitations as well, i.e. we have to know and fix the 3D shape of the disc beforehand.

The advantage of the firefly modelling is the lack of any grids. The flies are free to move within a predefined volume and their final distribution represents the location or distribution of emission within the disc. However, our approach does require that the disc has an axisymmetric homogenous density structure, as this is assumed for the computation of optical thickness effects that dominate the out-of-eclipse light curve shape.

### 5.2 Model Regularization

As in most of inversion problems, there are an infinite number of solutions if no regularization is employed. In order to overcome this, we have (like many previous authors) decided to search for the most axisymmetric image that can reproduce the observations (axisymmetric about the accretion disc rotation axis). In the case of the firefly approach this simply means a regularization term, where we have defined the disc in $`i`$ sectors, and counted the number of flies in each sector $`N_i`$ (here we have chosen to use 36 sectors). Then we maximize the entropy of such division i.e. we first compute

$`p_i={\displaystyle \frac{N_i}{_{i=1}^{36}N_i}}`$ (1)
where $`p_i`$is the fractional flux in each sector and the entropy term is

$`S={\displaystyle \underset{i=1}{\overset{36}{}}}p_i\mathrm{log}(p_i)`$ (2)
finally our merit function to be minimized is

$`F=\chi _\nu ^2\lambda S`$ (3)
where $`\chi _\nu ^2`$ is the standard reduced $`\chi ^2`$ term, $`\lambda `$ is the Lagrangian multiplier and $`S`$ the entropy term.

### 5.3 The fitting algorithm

We carry out the actual fitting using the same genetic algorithm variant as in Hakala et al. (2002). The only difference is that in this case the initial population of solutions for fitting consists of disc like swarms of random fireflies. These are then evolved in similar fashion to fit the eclipse profile and the overall light curve simultaneously. As the whole light curve has 532 points in B and 514 points in UV, we have decided to emphasize the eclipse phases (0.95-1.05) by summing the Chi-square values from this phase range three times. We typically use 50 individual swarms in each generation and after around 1500 generations, the fits stop evolving in a reasonable time.

In order to test our approach we have created a synthetic light curve using a known disc brightness distribution, added appropriate levels of noise and fitted that. The results of this exercise are shown in Figure 7. The top panel shows the whole orbital light curve together with the fit, the middle panel the phase range close to the eclipse phase and the bottom panel the synthetic brightness distribution (left) and the reconstructed distribution (right). The fits to the synthetic data are encouraging, with the reconstructed disc being close to the input disc. There is, however, a small amount of smearing, which is most likely caused by finite phase resolution, regularization and/or noise in the data.

### 5.4 The Results

We have plotted in Figures 8 and 9 the resulting brightness distribution from a fit to the B and UV band OM data in a case where no regularization was applied. The fit is excellent, but the resulting brightness distribution suffers slightly from the ‘scissors’ effects, i.e. the point spread function of the map has a scissors shape. This is because there is not enough information in the data to uniquely locate point like emission. Thus point source emission is spread out along the lines that mark the shadow of the secondary at the ingress and egress phases. This is a well known effect and was demonstrated by Horne (1985).

Using the regularization described earlier we obtain the fits for B and UV bands shown in Figures 10 and 11 respectively. The coordinate system for the fireflies is a Cartesian one, where the primary star is located in the origin and the X-Y plane is the orbital plane. The coordinates of the L1 point are roughly (0.7,0,0) and a ballistic freefall stream is overplotted. The observer lies ‘above’ the orbital plane, i.e. she/he would have a positive Z-coordinate.

Figure 12 shows the radial histograms of the distribution of flies. The resulting brightness distributions seem to indicate that whilst the UV emission comes mostly from the white dwarf or the boundary layer and the inner disc, the optical B band emission also contains a faint component from the accretion stream - disc interactions, i.e. the hot spot (and maybe even from the freefall part of the stream). This is also rather obvious directly from the light curve shapes, given the difference in eclipse assymmetry between the two passbands.

Furthermore, we find that the hot spot is located exactly on the ballistic trajectory, as expected. The distance is $``$ 0.3a from the white dwarf, which confirms the values reported by Wood et al. (1989). In addition to this we find marginal evidence for optical emission from the accretion stream prior to the main hot spot impact region. Our results show, that in the B band 90% of the emission comes within a distance of 0.29a from the white dwarf and 95% within 0.33a. This implies that the hot spot contributes roughly 5% of the total B band emission. The same radii in UV band are slightly greater in comparison (0.32a and 0.38a). At face value this implies that the B band radial emission distribution is concentrated more towards the white dwarf than in UV passband. This difference could be due to the fact that whilst our UV data was taken 4 days after an outburst, the B band data was taken more than a month later. Dwarf novae discs are thought to be more uniformly bright during the outbursts. Although the $`V`$ band mag was similar in both observations this does not preclude a difference in the size of the accretion disc in the two observations.

## 6 Discussion

### 6.1 The origin of the X-ray modulation

Ramsay et al (2001a) found that in X-rays the duration of the eclipse ingress and egress is 20–30sec, which is roughly the same time as that expected for the secondary star to eclipse the diameter of the white dwarf. These authors concluded that the X-ray emission region is hence extended in white dwarf longitude, probably centered near the spin equator. However, Wheatley & West (2003), using a re-analysis of the same data, conclude that the X-ray emission originated close to the upper pole on the white dwarf. In the latter scenario, an accretion column, or curtain, of material will rotate on the spin period of the white dwarf.

This is consistent by our spectral modelling in §4 which indicates that the most likely cause of these quasi-stable modulations is a variation in the absorption column over this period. If the 2193 sec period was the spin period of the white dwarf, how would that manifest itself as a variation in the absorption column? In the case of intermediate polars, some systems show an energy dependent X-ray variation at the spin period. This has been seen in EX Hya (Rosen, Mason & Cordova 1988, Cropper et al 2002) and is thought to be due to observing the emission site through an accretion curtain whose absorption column varies as the white dwarf rotates.

A second possibility may be due to a variation in the distance from the inner edge of the accretion disc to the white dwarf. Further, this distance could vary in azimuth in the binary frame. Thus the white dwarf-accretion disc interaction would lead to a variation at a beat period and could account for why the modulation is only seen for around 1/3 of the observation duration. This hypothesis requires that the 2193 sec period is the spin period of the white dwarf ($`\omega `$), whilst the 3510 sec modulation is the beat period between the spin and the 5436 sec orbital period ($`\mathrm{\Omega }`$). The beat period based on the 2193 s spin period would be ($`\omega \mathrm{\Omega }`$) 3676 sec, which is 6.6$`\sigma `$ distant from the best fit period. However, if the spin period is actually 1.7$`\sigma `$ shorter than the best fit value (giving 2133 sec) then the observed 3510 sec period would be exactly the expected $`\omega \mathrm{\Omega }`$ beat period. In this scenario, any very short periods (for instance the 17.94 and 18.15 sec periods reported by Marsh & Horne, 1998) would have to be explained in other terms, for instance white dwarf oscillations.

Finally, Marsh & Horne (1998) suggested that one of the $``$18 sec periods could represent $`\omega `$, or alternatively that one of the periods was the Keplerian period and the other was the beat period between this and the spin period (implying a spin period of $``$1550 sec). If $``$1550 sec was the spin period, then most of the amplitude peaks in our Figure 1 could be explained as naturally occuring sidebands of the spin and orbital frequencies in the following manner: 1248 sec ($`\omega +\mathrm{\Omega }`$), 1574 sec ($`\omega `$), 2193 sec ($`\omega \mathrm{\Omega }`$), 4844 sec (2($`\omega \mathrm{\Omega }`$)). The main drawback with this hypothesis is that there is no simple explanation for the very significant 3510 sec period. Also, it is possible that not all of the 5 periods mentioned above are in fact physical (cf §3).

Recently, the short period variability in cataclysmic variables has been extensively discussed in a series of papers by Warner & Woudt (most recently in Warner, Woudt & Pretorius (2003)). In conclusion they have found that typically $`P_{QPO}15P_{DNO}`$. In addition they also find second type of DNO’s with: $`P_{lpDNO}4P_{DNO}`$. Assuming that one of the periods seen here, 2193s or 3500s, would be a QPO period the above mentioned relationship would then predict DNO periods around 150s and 230s respectively. There is no evidence for such periods in our data. However, this is expected since the DNO and QPO periods are normally found during outbursts, and during our observation the source was in quiescence.

We conclude that based on our data, we cannot provide a unique explanation for the 2193 sec and 3510 sec periods. We can, however, conclude that these periods occur due to changes in the line of sight absorbing column.

### 6.2 The Optical Light Curve Modelling

Much previous work has been done on the eclipse profiles of OY Car ranging from infrared to X-rays (see for instance Berriman 1984, Cook 1985, Krzeminski & Vogt 1985, Schoembs 1986, Berriman 1987, Wood et al. 1989, Rutten et al. 1992, Horne et al 1994, Bruch, Beele and Baptista 1996, Pratt et al. 1999, Ramsay et al. 2001a, Wheatley & West 2003). Wood et al. (1989) use the optical eclipse profiles to measure accurate system parameters and assuming a ballistic free fall trajectory for the accretion stream they deduce a disc outer radius of approximately 0.31a. Rutten et al. (1992) re-analysed the set of eclipse profiles from Vogt (1983) that cover the rise into an outburst. Their eclipse mapping analysis shows that whilst the quiescence eclipse profiles are dominated by the white dwarf and the hot spot, the rest of the disc becomes dominant in the outburst. The HST UV observations by Horne et al. (1994) demonstrate that the eclipse profiles become more symmetric and ‘tophat shaped’ towards the shorter wavelengths. This indicates that the white dwarf and/or the boundary layer is the major source of emission in the UV, whilst the hot spot contribution rises towards the optical wavelengths.

Our light curve modelling confirms the previous results on the origin of the quiescent eclipse profiles. In particular they indicate, that both in B and UV bands, the majority of the emission arises close to the white dwarf, possibly either the boundary layer or the white dwarf itself. This is also supported by the ingress and egress time analysis of the the same B, UV dataset by Ramsay et al. (2001a).

The freefall accretion stream and especially the hot spot, where the stream hits the disc outer edge is mainly emitting in the B band and hardly visible in the UV maps. Our modelling also suggests that in addition to the hot spot, the ballistic accretion stream might contribute to the B band emission as well.

Our novel approach allows simultaneous fitting of the eclipse profiles and orbital light curves of cataclysmic variables. This removes bias introduced by artificially having to remove the orbital modulation prior to the eclipse profile mapping. Furthermore, the ‘firefly approach’ has potential to be easily converted into full 3D mapping, if a proper regularisation can be employed.

## 7 Conclusions

We have reanalysed the XMM-Newton observations of the dwarf novae OY Car. In particular, we have reanalysed the X-ray light curve and spectrum and developed a new model to fit the OM B band and UV data.

We have investigated the X-ray light curve over the time interval where the previously reported 2240 sec modulation was strongly seen. This has allowed us to refine this period to 2193 sec. We also detect a modulation at a period of 3510 sec which we also believe to be significant. We have performed phase resolved X-ray spectroscopy on these periods and find that the soft X-ray light curve is anti-correlated with the absorption column density. This implies that the light curve modulation is due to us observing the X-ray emission region through a rotating accretion column or curtain. Although, not conclusive, this supports the hypothesis of Ramsay et al (2001a) that the 2193 sec is the spin period of the white dwarf. Furthermore, this is in agreement with findings of Wheatley & West (2003), where they conclude that the X-ray emission probably originates near the polar region on the white dwarf surface.

To model the OM B and UV light curves, we have applied the fireflies model of Hakala et al (2002). We simulated light curves using a given accretion disc brightness distribution and used our technique to fit the synthetic light curves. The resulting distribution of fireflies is encouraging. We then applied the model to the OM data and find that that the hot-spot, where the accretion stream hits the disc, is exactly on the ballistic trajectory. Further, we find some evidence that the B band emission is more concentrated than the UV emission, although this maybe due to the fact that the B band data were taken 1 month after the UV data. Our technique can be applied to other disc accreting systems.

###### Acknowledgements.

This paper is based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and the USA (NASA). PJH is an Academy of Finland research fellow.


format: text

text:
# Acknowledgments

## Acknowledgments

It is a pleasure to thank Gabriele Ferretti and Rainer Heise for a fruitful collaboration on this and related work, and the organizers of the Johns Hopkins Workshop in Göteborg and the RTN Workshop in Copenhagen for giving me the opportunity to present it, and for the nice atmosphere. This work was and is partly supported by EU contracts HPRN-CT-2000-00122 and HPRN-CT-00131, by the “Actions de Recherche Concertées” of the “Direction de la Recherche Scientifique - Communauté Française de Belgique”, by a “Pôle d’Attraction Interuniversitaire” (Belgium) and by IISN-Belgium (convention 4.4505.86). R.A. is a Postdoctoral Researcher of the Fonds National de la Recherche Scientifique (Belgium).


text:
# Inherent size constraints on prokaryote gene networks due to “accelerating” growth

## I Introduction

The rapidly expanding field of network analysis, reviewed in Dorogovtsev\_02\_10 ; Albert\_02\_47 , has provided examples of networks exhibiting “accelerating” network growth, where link number grows faster than linearly with network size. For instance, the Internet Faloutsos\_99\_25  appears to grow by adding links more quickly than sites though the relative change over time is small and the Internet appears to remain scale free and well characterized by stationary statistics Vasquez\_02\_066130 . Similarly, the number of links per substrate in the metabolic networks of organisms appears to increase linearly with substrate number Jeong\_00\_65 , the average number of links per scientist in collaboration networks increases linearly over time Dorogovtsev\_00\_33 ; Vasquez\_00\_0006132 ; Barabasi\_01\_0104162 ; Barabasi\_02\_590 ; Vasquez\_03\_056104 , and languages appear to evolve via accelerated growth Dorogovtsev\_01\_2603 .

In the main, the chief focus of these studies has been on locating parameter regimes allowing accelerating networks to maintain scale free statistics and thereby to allow continued unconstrained growth. For example, an early study considered a growing network receiving $`N^\alpha `$ new links for $`\alpha >0`$ when the network size is at $`N`$ nodes, but restricted analysis to the case $`\alpha 1`$ as “Obviously, $`\alpha `$ cannot exceed 1 (the total number of links has to be smaller than $`N^2/2`$ since one may forbid multiple links).” and “The density of connections in real networks remains rather low all the time, so one may reasonably assume that $`\alpha `$ is small.” Dorogovtsev\_01\_025101 . Equivalent limits were considered in Ref. Sen\_0310513 . In such restricted parameter regimes networks could maintain scale free statistics, though this result carries the implicit but unexamined finding that alternate parameter regimes permit transitions from stationary to nonstationary statistics. This paper builds on these implicit findings.

Accelerating networks are more prevalent and important in society and in biology than is commonly realized—see the survey in Ref. Gagen\_03\_accel\_survey . In fact, any network that requires functional integration and organization (where the activity of any given node is dependent on the state of the network or different subnetworks) is by definition an accelerating network, that is, as the network expands, the proportion of the network devoted to control and regulation expands disproportionately. This in turn means that all such networks, sooner or later, must be limited in their size and complexity, which limitations can only be breached by changing either the physical nature of the control architecture (a state transition) or by reducing the functional integration. In the latter case, where networks are hitting a complexity limit, further growth in network size will likely display structural transitions from randomly connected, to scale free statistics, to densely connected and perhaps finally to fully connected statistics. Should such networks be unable to successfully complete these transitions for any reason, then it is likely that network growth must cease entirely or that either a transition to a nonaccelerating structure is required to permit further growth or novel technologies must appear allowing the continuation of accelerated growth. Exemplar accelerating networks displaying such size limits or structural transitions include (a) all forms of economic markets where the latest price offered by any participant instantly affects all other participants, (b) industrial companies and sectors implementing a Just-In-Time business model where any worker can halt the entire production system, (c) error propagation networks linking an error source with all affected nodes as studied in software analysis and in models of the propagation of diseases, bushfires, cracks, and electricity grid failures, (d) in any dynamical system dependent on relative quantities so changes in one node instantly affects every other node such as relative transcription factor binding probabilities or relative evolutionary fitness, (e) in computer hardware and in cluster and grid supercomputer networks, and (f) in organizational networks Gagen\_03\_accel\_survey . In fact, it is well understood that social networks only take on small world statistics when the network is large enough—in small towns everyone one knows everyone else so social networks are accelerating, and social networks make a transition to small world statistics only as individual nodes saturate their connectivity limits Watts\_99\_493 . Similar observations can be made about the scale free Internet and World Wide Web—when sufficiently small, these networks were likely accelerating until connectivity capacities saturated forcing a transition to scale free structures to permit further growth Gagen\_03\_accel\_survey .

This paper develops an accelerating network model of prokaryotic (single celled) gene regulatory networks to investigate size and complexity limits inherent in the adoption of an accelerating architecture. Because our focus is on structural transitions, we explicitly do not need to restrict the degree of acceleration to low values of $`\alpha 0`$. Rather, we permit this parameter to take on any value including $`\alpha >1`$ and ensure that the network is not saturated by making link formation probabilistic. The resulting novel “probabilistic” accelerating networks grow by adding on average $`pN^\alpha `$ new links with $`\alpha >0`$ and otherwise arbitrary provided the probability of adding a link is suitably constrained $`p1`$ so that total link number remains less than of order $`N^2`$.

The gene regulatory model presented here is motivated by comparative genomics findings that the total number of regulatory proteins controlling gene expression (links) scales quadratically with the number of genes or operons (nodes) in prokaryotes vanNimwegen\_03\_479 ; Croft\_03\_unpub . This quadratic growth results as the number of links made by a regulator exploiting homology dependent (sequence specific) interactions scales proportionally to the number of randomly drifting promotor sequences or effectively, with gene number Croft\_03\_unpub . Hence, gene regulatory networks are inherently accelerating—the probable number of links per regulator $`pN^\alpha `$ increases linearly with node number with $`\alpha =1`$, so consequently, the total number of links scales quadratically as $`pN^{\alpha +1}`$. In small and sparsely connected networks, most links come from different regulators suggesting that regulator number also scales quadratically with gene number, $`pN^{\alpha +1}`$. Such an accelerating network would be characterized initially by sparse connectivity at low gene numbers and subsequently by denser connectivity at high gene numbers as networks attempt a transition to a densely connected regime. If the evolving networks can successfully make this transition, the evolutionary record will display a transition in network statistics for some critical network size $`N_c`$. Conversely, if these networks, optimized by evolution in the sparse regime, are unable to make the transition to the densely connected regime, the evolutionary record would show a strict size limit $`NN_c`$ at some critical network size. But this is exactly what is observed. All prokaryotic gene numbers and genomes are indeed of restricted size (less than about 10,000 genes with genomes of between 0.5 and 10 megabases Casjens\_98\_33 ), in contrast to the genomes of multicellular eukaryotes (with for humans, about 30,000 genes and a genome of about 3 gigabases Int\_Human\_genome\_01\_86 ; Venter\_01\_13 ). Ref. Croft\_03\_unpub  predicted the size limit $`N_c20,000`$ genes as continued genome growth requires the number of new regulators to exceed the number of nonregulatory nodes.

A satisfactory model of prokaryotic gene regulatory networks requires some novel features. As mentioned above, we introduce probabilistic link formation to allow rapid accelerated growth and correspondingly stricter size limits. (A different but related mechanism was introduced in Refs. Liu\_02\_036112 ; Goh\_02\_108701  which considered the effects of stochastic fluctuations in the number of added links with each additional node.) In addition, we employ directed links and partition nodes into two classes where “regulators” can source outbound regulatory links to regulate other nodes (both regulators and non-regulators), while “non-regulators” cannot source outbound links. (Ref. Cheng\_02\_066115  has previously considered networks of distinguishable nodes.) Further, experimental evidence presented below indicates that the distribution of inbound links is compact and exponential while the distribution of outbound links is long-tailed and likely scale-free. As a result, the heads and tails of our directed links are placed according to two distinct distributions. Altogether, these features allow the reproduction of the observed features of prokaryote gene regulatory networks and satisfactorily predicts the maximum prokaryotic gene count.

Our approach reproducing accelerating network statistics for growing prokaryote genomes complements and informs alternate networking approaches seeking to deduce or simulate the regulatory networks of particular organisms from gene perturbation and microarray experiments Wagner\_01\_1183 ; Bhan\_02\_1486 ; Lukashin\_03\_1909 .

In Section II we canvass the available literature to characterize the statistics of prokaryote gene regulatory networks. This then allows the construction of accelerating growth network models in Section III where we use the continuous approximation and simulations to analyze network statistics. The size constraints inherent in accelerating prokaryote regulatory networks are modelled in Section IV.

## II Overview of prokaryote gene networks

Ongoing genome projects are now providing sufficient data to usefully constrain analysis of the gene regulatory networks of the simpler organisms. Ref. vanNimwegen\_03\_479  first noted the essentially quadratic growth in the class of transcriptional regulators ($`R`$) with the number of genes ($`N_g`$) in bacteria with the observed results

$$R\{\begin{array}{cc}N_g^{1.87\pm 0.13},\hfill & \text{transcriptional regulation}\hfill \\ & \\ N_g^{2.07\pm 0.21},\hfill & \text{two component systems}\hfill \\ & \\ N_g^{2.03\pm 0.13},\hfill & \text{transcriptional regulation}\hfill \\ & \\ N_g^{2.16\pm 0.26},\hfill & \text{transcriptional regulation}.\hfill \end{array}$$
(1)
Here, the top two lines refer to different classes of regulators while the bottom two lines are the results of a crosschecking analysis of two alternate databases. Quoted intervals reflect 99% confidence limits vanNimwegen\_03\_479 . The explanation for this quadratic growth was that each additional transcription factor doubles the number of available dynamical states which, it was posited, allows for a doubling in the fixation probabilities for this class of genes.

As noted above, Ref. Croft\_03\_unpub  provides an alternate theoretical analysis predicting quadratic growth in any regulatory network exploiting homology dependent interactions and analyzed 89 bacterial and archeael genomes to determine the relations

$$R=\{\begin{array}{cc}aN_g^b=(1.6\pm 0.8)10^5N_g^{1.96\pm 0.15}\hfill & (r^2=0.88)\hfill \\ & \\ pN_g^2=(1.10\pm 0.06)10^5N_g^2\hfill & (r^2=0.87)\hfill \\ & \\ cN_g=(0.055\pm 0.004)N_g\hfill & (r^2=0.75).\hfill \end{array}$$
(2)
In this paper, accelerating networks will be based on the quadratic second line (while nonaccelerating models presented in later work will work with the linear third line Gagen\_0312022 ). In all cases, the limits reflect 95% confidence levels. For completeness, the data is shown in Fig. 1. The observed quadratic growth implies an ever growing regulatory overhead so there will eventually come a point where continued genome growth requires the number of new regulators to exceed the number of nonregulatory nodes, and based on this, Ref. Croft\_03\_unpub  predicted an upper size limit of about 20,000 genes, within a factor of two of the observed ceiling.

Earlier surveys of bacterial genomes noted that larger genomes harboured more transcription factors per gene than smaller ones Cases\_03\_248 , with this trend attributed to the need in larger genomes for a more complex network of regulatory proteins to achieve coordinated expression of a larger set of cellular functions, and to selection in complex environments leading to enrichment in transcription factors allowing regulation of gene expression and signal integration. A similar upward trend in the proportion of regulators as a fraction of genome size with increasing genome size was observed in Ref. Stover\_00\_959  attributed to a need for an increasing responsiveness in diverse environments, with confirming observations in Ref. Bentley\_02\_141 .

Prokaryotes typically group their DNA encoded genes in operons, co-regulated functional modules of average size 1.70 genes each in E. coli which value we treat as typical though in reality, operon size decreases slightly with genome size Cherry\_03\_40 . Each operon can be either unregulated and so constitutively or stochastically expressed or subject to combinatoric regulation by multiple regulatory protein transcription factors binding to each operon’s promotor sequence.

Again assuming that E. coli is typical, any given regulatory protein affects an average of about 5 operons with this distribution being long tailed Shen\_Orr\_02\_64  so the majority of regulators affect only one operon though some regulators (CRP) can affect up to 71 operons or 133 genes Thieffry\_98\_43 . (This latter reference estimated that each regulator controls on average 3 genes.) More recent estimates have the transcription factor CRP, a global sensor of food levels in the environment, regulating up to 197 genes directly and a further 113 genes indirectly via 18 other transcription factors MandanBabu\_03\_1234 . (To observe the long tailed distribution, see Fig. 2 of Ref. Thieffry\_98\_43  and Fig. 4 of Ref. MandanBabu\_03\_1234 .)

However, the number of inputs taken by an operon is characterized by a compact exponential distribution with a rapidly decaying tail so the majority of regulated operons are controlled by a single regulator while very few regulated operons are controlled by four, five, six or seven regulators Thieffry\_98\_43 ; Shen\_Orr\_02\_64 ; MandanBabu\_03\_1234 . In particular, Ref. Thieffry\_98\_43  examined 500 regulatory links from about 100 regulators to almost 300 operons to estimate that each regulated operon takes on average 2 inputs though Fig. 2 of this reference suggests an average input number of about 1.5. Similarly, Ref. Shen\_Orr\_02\_64  suggests that 424 regulated operons receive 577 links giving an average input number of 1.4, while Ref. MandanBabu\_03\_1234  estimates that 327 regulated operons receive 524 links giving an average input number of 1.6.

## III Accelerating prokaryote network models

We extend the gene network model of Ref. Thieffry\_98\_43  to construct an accelerating network model of prokaryote regulatory gene networks. Prokaryotes typically pack their $`N_g`$ genes into a lesser number of $`N=N_g/g_o`$ co-regulated operons where we assume that operons contain exactly $`g_o=1.70`$ genes. Of the existing operons, $`O_r`$ are regulated operons and $`O_u=NO_r`$ are unregulated operons. Of the total number of operons, there are $`R`$ regulatory operons whose regulatory interactions are directed links from regulatory operons to regulated operons. Under the assumption that there is only one regulatory gene per regulatory operon, the observed quadratic relation of Eq. 2 becomes

$$R=pg_o^2N^2.$$
(3)
When regulators and regulatory links are very rare, i.e. when genomes are small, it is likely that every new link is associated with a new regulator so the number of links varies roughly quadratically with operon number. We write

$$L=lN^2,$$
(4)
where $`l`$ denotes the probability of forming a particular beneficial link per operon. The value for $`l`$ will be approximately $`pg_o^2`$, but the exact relation must be derived from the details of the implemented model.

Each regulatory link between nodes is directed, and characterized by two distinct distributions describing respectively the placement of the heads and tails of each link. Only a relatively few nodes are regulatory, and of these, the number of outbound link tails per regulatory node are described by a size dependent long-tailed distribution with average about $`t5`$. Such a long-tailed distribution requires that link tails be preferentially attached to an existing regulatory operon or equivalently, the associated regulated operon must possess one promotor binding site (among others) that binds that particular regulator. Consequently, the preferential selection of regulators means that the promotor sequences of newly regulated nodes cannot be randomly chosen—randomly drifting promotor sequences would be as likely to match any one regulator as another. A plausible physical explanation for the preferential attachment of link tails to existing regulators is that newly fixated operons come largely from gene duplication events Yanai\_00\_2641  where some of the duplicated promotor binding sites are under strong selective constraint while other binding sites and the operon genes can drift freely. Gene duplication then implies that in a genome of size $`N`$ operons, if some regulator $`n_j`$ has $`t_{jN}`$ outbound regulatory links to approximately $`t_{jN}`$ regulated operons, then the probability that a newly fixated operon is also regulated by $`n_j`$ is simply the proportion of such regulated operons in the genome, or $`t_{jN}/N`$. This implements the required preferential attachment as the resulting rate of growth in the number of links attached to node $`n_j`$ is also then proportional to $`t_{jN}`$. If there is also some probability of the appearance of novel promotor sequences, these combined processes suffice to produce the observed scale free distributions. This model is roughly consistent with recent estimates of the relative contributions to prokaryote genome growth which suggest that horizontal gene transfer rates $`\gamma _h`$ are roughly one third of gene loss rates $`\gamma _h=\gamma _l/3`$ and roughly one half of vertical inheritance or gene genesis rates $`\gamma _h=\gamma _v/2`$ leading to roughly constant sized genomes over long times (as $`\dot{N}\gamma _h+\gamma _v\gamma _l0`$), while “it is remarkable that phylogenetic distributions of at least 60% \[and up to 75%\] of protein families can be explained merely by vertical inheritance.” Kunin\_03\_1589 . Similarly, three quarters of examined transcription factors in Ref. MandanBabu\_03\_1234  were two-domain proteins with shared domain architectures leading to the estimate that about 75% of transcription factors have arisen as a consequence of gene duplication (though the joint duplication of regulatory regions and of regulated genes or of transcription factors together with regulated genes is more rare). A further implication of these gene duplication processes is that, in the main, regulators can only appear on entry to the genome—a potential regulator lacking any target matches in a given genome will never form any links when most operons arise from promotor preserving duplication events. This allows us to considerably simplify our model, and hereafter, we only allow regulators to appear on their entry to the genome. Of course, more realistic but considerably more complicated models are possible.

In contrast to the relatively small number of regulatory nodes, all nodes can themselves be regulated by inbound links and in fact, can be multiply regulated as promotor regions can contain more than one binding site. Further, the many used and unused promotor region binding sites broadly sample the space of possible binding sites so only a small fraction of nodes will be regulated by any one regulator. As a result, the number of inbound link heads per node is described by a size dependent exponential distribution with a low average of $`h1.5`$ as typically results from the random or non-preferential attachment of inbound links to operon promotor sequences.

We suppose that the operon network grows by the sequential addition of numbered nodes $`n_k`$ for $`1kN`$, and that at network size $`k`$, node $`n_i`$ ($`1ik`$) has $`t_{ik}`$ outbound tails and $`h_{ik}`$ inbound heads. We do not model the many trials of potential genes over many generations and merely include fixated genes in our count—that is, drifting sequence is not counted as part of the fixated genome. This further implies the sequence of established nodes is under severe selective constraint and unable to drift so consequently new links cannot be added between existing nodes. (If a proportion $`fN`$ of existing nodes can explore novel sequence space in time $`dt`$, then the number of new regulators increases as $`dRfN^2dt`$, and as $`N`$ is itself a function of time, this integrates to generate a non-quadratic relation between regulator and operon number which is not observed.)

For clarity, Fig. 2 preempts later calculations and depicts a statistically generated version of an E. coli genome where nodes are placed sequentially counterclockwise in a circle (for convenience only). Alternative genome models may be distinguished by the age distribution of regulators, regulated operons and their link numbers, and these are indicated in this figure. In particular, Fig. 2 shows a highly nonuniform distribution of regulators and outbound link numbers with gene age in contrast to a uniform distribution of regulated operons and of inbound link numbers. (It will turn out that these latter age-independent distribution are only present when regulator number grows quadratically with genome size.)

These distributions result from the physical processes underlying the formation of regulatory links in prokaryotes. As discussed above, a substantial proportion of the gene regulation network of prokaryotes is enacted via homology dependent interactions as when sequence specified protein transcription factors bind to specific promoter sequences. The undirected nature of evolutionary searches means that gene regulatory networks fundamentally exploit the same sequence matching algorithms used in comparative genetics where the probability of obtaining matches between a single given trial sequence of some small fixed length and an entire genome scales proportionately to genome length—doubling genome length doubles the probability of a match. Hence, the expected number of links formed per regulator scales linearly with present genome size. As the number of source trial sequences also scales with genome length, the expected number of matches between all regulators and all regulated operons scales quadratically with genome length, or effectively, with operon number assuming constant sized operons over the evolutionary record.

As a consequence, on entry into the genome, each new gene has some probability of being a regulator dependent firstly on its suitability to bind DNA and secondly on the linearly increasing expected number of acceptable binding targets present in the genome on entry (or at later times). As discussed above, the predominance of vertical gene genesis events allows a simplified model wherein the probability of a new node being regulatory is determined solely by the number of available links present at the time of entry. We assume then that on entry into the genome each new node $`n_k`$ can form $`2k1`$ links with nodes $`n_1,\mathrm{},n_k`$ consisting of a single self-regulatory link from node $`n_k`$ to itself with probability $`l`$, $`(k1)`$ regulatory outbound links to the existing nodes each with equal probability $`l`$, and, provided that sufficient regulators already exist, $`l(k1)`$ inbound regulatory links from some subset of the existing regulators chosen according to preferential attachment. (For consistency, we can only add $`lk`$ distinct regulatory links to node $`n_k`$ provided there are at least this many regulators in existence. From Eq. 3, the average number of regulators $`pg_o^2k^2`$ must be greater than the number of regulatory links $`lk`$, and this will be satisfied for $`k>l/(pg_o^2)1`$.) As a result, the total number of heads or tails attached to node $`n_k`$ on entry to the genome ranges between $`0`$ and $`k`$, with each link formed with probability $`l`$. Hence, the respective probabilities that the initial number of heads $`h_{kk}=j`$ or the initial number of tails $`t_{kk}=j`$ for node $`n_k`$ is

$$P(j,k)=\left(\genfrac{}{}{0pt}{}{k}{j}\right)l^j(1l)^{kj},$$
(5)
with the proviso that all the inbound links can only be added to node $`n_k`$ if there is a sufficient number of regulators among the nodes $`n_1,\mathrm{},n_k`$. The average number of inbound and outbound links is identical, $`t_{kk}=h_{kk}=lk`$ showing linear growth in link number with increasing network size. The addition of node $`n_k`$ and its links will increase the probable number of heads attached to earlier nodes $`n_j`$ for $`1j(k1)`$ so $`h_{jk}h_{jj}`$, while the probable number of tails outbound from node $`n_j`$ increases $`t_{jk}t_{jj}`$ if and only if that node is regulatory with $`t_{jj}>0`$.

As regulators can only be created on entry to the genome, the distribution of regulators at any time is specified by the distribution $`P(j,k)`$ for $`t_{kk}`$. Using Eq. 5, the probability that node $`n_k`$ is a regulator is $`1P(0,k)`$, so for a network of $`N`$ nodes, the predicted total number of regulators is

$`R`$ $`=`$ $`{\displaystyle \underset{k=1}{\overset{N}{}}}\left[1(1l)^k\right]`$ (6)
$`=`$ $`N{\displaystyle \frac{1l(1l)^{N+1}}{l}}`$
$``$ $`{\displaystyle \frac{l}{2}}N(N+1).`$
The exact top line shows the expected behaviour for the number of regulators in the respective limits $`l0`$ giving $`R0`$, and $`l1`$ giving $`RN`$. The approximate relation in the third line can be compared to the observed Eq. 3 and immediately suggests $`l2pg_o^2`$, while a fit to the more accurate top line gives the connection probability as

$$l=1.15\times 2pg_o^2=\mathrm{\hspace{0.33em}7.31}\times 10^5.$$
(7)
This probability value suggests an average promotor binding site length of $`\mathrm{log}_4l=6.9`$ bases. The average number of links per regulator using the second line of Eq. 6 is then approximately $`L/R2`$, while the more accurate top line with $`N=2528`$ operons for E. coli Cherry\_03\_40  gives $`L/R=2.12`$, about a factor of two from the observed value of $`5`$ for E. coli Shen\_Orr\_02\_64 .

### III.1 Random distribution of regulated operons

The distribution of link heads for all nodes (with possession of a link head designating a regulated node), can be straightforwardly calculated under the assumption that the $`t_{kk}lk`$ new tails added with node $`n_k`$ are randomly distributed across the $`k`$ existing nodes so on average, each existing node receives $`l`$ links. To build insight, it is useful to consider the general case where $`t_{kk}h_{kk}lk^\alpha `$ for $`\alpha 0`$. Setting $`\alpha =0`$ adds with some probability a constant number of links with each new node, $`\alpha =1`$ adds a linearly growing number of probable links with each new node, $`\alpha =2`$ adds a quadratically growing number of probable links with each new node, and so on. The total number of links present in the network is then

$$_0^N2lk^\alpha =\frac{2lN^{\alpha +1}}{\alpha +1}$$
(8)
The continuous approximation Barabasi\_99\_17 ; Barabasi\_99\_50 ; Dorogovtsev\_01\_056125  for links randomly distributed over $`k`$ existing nodes determines the number of inbound head links for node $`n_j`$ according to

$`{\displaystyle \frac{h_{jk}}{k}}`$ $`=`$ $`{\displaystyle \frac{t_{kk}}{k}}`$ (9)
$`=`$ $`lk^{\alpha 1}.`$
This can be integrated with initial conditions $`h_{jj}lj^\alpha `$ at time $`j`$ and final conditions $`t_{jN}lN^\alpha `$ at time $`N`$ to give

$$h_{jN}=\{\begin{array}{cc}l+l\mathrm{ln}\frac{N}{j}& \text{if }\alpha =0\\ & \\ \frac{l}{\alpha }N^\alpha +\frac{l(\alpha 1)}{\alpha }j^\alpha & \text{if }\alpha >0.\end{array}$$
(10)
Integration of these link numbers over all node numbers $`j`$ gives the required total number of links as in Eq. 8. For $`0\alpha <1`$, the number of links per node is monotonically decreasing with node number. However, for $`\alpha =1`$ and only in this case, the final distribution is independent of node number $`j`$ because earlier nodes receive exactly enough links from latter nodes to balance the initially biased distribution of heads $`h_{jj}lj`$, so in the end, all nodes receive on average the same number of inbound regulatory links $`h_{jN}=lN`$ for $`1jN`$. For faster acceleration rates, $`\alpha >1`$, the number of links per node is monotonically increasing as later nodes receive a greater number of links on entry to the genome and this imbalance is not corrected.

The possibility of monotonically increasing numbers of links with node number in accelerating networks has not previously been considered. This possibility requires modifying the usual continuum approach Barabasi\_99\_17 ; Barabasi\_99\_50 ; Dorogovtsev\_01\_056125  so the final inbound link distribution is obtained via

$`H(k,N)`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle _0^N}𝑑j\delta (kh_{jN})`$ (11)
$`=`$ $`\pm {\displaystyle \frac{1}{N}}\left({\displaystyle \frac{h_{jN}}{j}}\right)^1\text{at }[j=j(k,N)],`$
where $`j(k,N)`$ is the solution of the equation $`k=h_{jN}`$. The top line is used when all nodes possess the same average link number while the second line is applicable with the plus (negative) sign when the average numbers of links per node is monotonically increasing (decreasing) with node number. Non-monotonic cases require alternate approaches.

Under quadratic growth in total link number when $`\alpha =1`$, and only in this case, the final distribution of link heads is independent of node number and evaluated using Eq. 11 to give

$`H(k,N)`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle _0^N}𝑑j\delta (klN)`$ (12)
$`=`$ $`\delta (klN).`$
As expected, a compact final link distribution results when all nodes have an average of $`t_{jN}=lN`$ inbound regulatory links at time $`N`$. This distribution calculated under the continuous approximation equates to one where in reality, each node receives a controlling head with probability $`l`$ from every other node (though in practise, the total number of received links is of order unity). Hence, for any node in a network of size $`N`$, the actual probability of having $`k`$ heads is

$$H(k,N)=\left(\genfrac{}{}{0pt}{}{N}{k}\right)l^k(1l)^{Nk}.$$
(13)
A network simulation with linear growth of link numbers per node model (Eq. 5) serves to validate this predicted final distribution. Fig. 3 compares the predicted distribution $`H(k,N)`$ against observed distributions for typical simulated networks of various sizes with negligible discrepancies.

For a network of size $`N`$, the probability that any given operon is unregulated is $`H(0,N)`$ so the expected number of unregulated operons summed over all $`N`$ nodes is

$$O_u=N(1l)^N.$$
(14)
This determines the number of regulated operons as

$$O_r=NO_u=N\left[1(1l)^N\right],$$
(15)
showing the expected behaviour as $`l1`$ giving $`O_rN`$ and $`l0`$ giving $`O_rlN^2=L`$ as each of the sparsely distributed links hits a distinct regulated operon. We note that random gene duplication and deletion events will not change the $`H(k,N)`$ distribution (other than changing $`N`$) as all nodes are identically connected on average. The $`H(k,N)`$ distribution appears in Fig. 2 which shows a uniform (age-independent) distribution of regulated nodes over the genome, and this uniformity is only expected for $`\alpha =1`$ corresponding to linear growth in link numbers per node and quadratic growth in regulator numbers.

These predictions can be compared to observation. For the E. coli network of size $`N=2528`$ operons or 4289 genes Cherry\_03\_40 , the predicted proportion of regulated operons receiving $`k>0`$ inputs is

$$P_h(k)=\frac{H(k,N)}{1H(0,N)},$$
(16)
and is shown in Fig. 4. Here, the calculated distribution closely approximates the compact exponential distribution observed for E. coli shown in Fig. 2(d) of Ref. Cherry\_03\_40  and of Fig. 5 of Ref. MandanBabu\_03\_1234 , though it underestimates the numbers of regulated operons with 4, 5, 6 and 7 inputs—essentially no regulators are predicted to have 5 or more inputs for genomes of size $`N=2528`$ operons. In addition, the average number of inbound regulatory links per operon (for all operons) is $`k=L/N=lN=0.19`$, while the average number of inbound regulatory links for regulated operons is $`k_r=L/O_r1`$. A more accurate calculation using the specific values for E. coli gives $`k_r=L/O_r=1.10`$, very close to the E. coli value of 1.5 or 1.6 noted in Refs. Thieffry\_98\_43 ; Shen\_Orr\_02\_64 ; MandanBabu\_03\_1234 .

### III.2 Scale-free distribution of regulator operons

On entry into the genome, node $`n_k`$ sources on average $`lk`$ outbound regulatory links and this linear growth in link number means that more recent nodes are more likely to be immediately regulatory and more likely to be highly connected on genome entry. However, node $`n_k`$ will also receive on average $`lk`$ inbound regulatory links whose tails will be preferential attached to existing regulators. The final distribution of link number with age will depend on the rate at which earlier nodes under preferential attachment can attract links relative to the linearly increasing link numbers of later regulatory nodes.

On entry at time $`k`$, node $`n_k`$ receives $`h_{kk}lk`$ inbound links from existing regulatory nodes in the set $`n_1,\mathrm{},n_k`$. As previously, we gain insight by considering the general case where $`t_{kk}h_{kk}lk^\alpha `$ for $`\alpha 0`$ (though we continue to use the distribution $`P(j,k)`$ of Eq. 5 to determine both the number of links $`j`$ prior to exponentiation and regulatory probability so consequently the number of regulators continues to increase quadratically according to Eq. 6). As a result, the need to ensure that all regulatory links to node $`n_k`$ are distinct requires that new link number $`lk^\alpha `$ be less than the number of existing regulators $`lk^2`$ requires $`\alpha 2`$. The $`h_{kk}`$ new tails added with node $`n_k`$ are preferentially attached to the existing regulatory nodes $`n_j`$ with probability proportional to the number of existing regulatory links for that node at time $`k`$, i.e. $`t_{jk}`$. Using the continuous approximation Barabasi\_99\_17 ; Barabasi\_99\_50 ; Dorogovtsev\_01\_056125 , the rate of growth in outbound link number for node $`n_j`$ is then approximately

$$\frac{t_{jk}}{k}=h_{kk}\frac{t_{jk}}{_0^kt_{jk}𝑑j}.$$
(17)
The denominator here is a probability weighting to ensure normalization and is the total number of outbound links for all nodes. Following Dorogovtsev\_02\_10 , we can evaluate the denominator using the identity

$$\frac{}{k}_0^kt_{jk}𝑑j=_0^k\frac{}{k}t_{jk}𝑑j+t_{kk}.$$
(18)
This can be evaluated using Eq. 17 noting $`t_{kk}h_{kk}lk^\alpha `$ giving

$$\frac{}{k}_0^kt_{jk}𝑑j=2lk^\alpha ,$$
(19)
which can be integrated determining the denominator of Eq. 17 to be

$$_0^kt_{jk}𝑑j=\frac{2l}{\alpha +1}k^{\alpha +1}.$$
(20)
This is in agreement with Eq. 8. Substituting this value into Eq. 17 gives

$$\frac{t_{jk}}{k}=\frac{\alpha +1}{2}\frac{t_{jk}}{k}.$$
(21)
Finally, this can be integrated with initial conditions $`t_{jj}lj^\alpha `$ at time $`j`$ and final conditions $`t_{jN}`$ at time $`N`$ to give

$$t_{jN}=lN^{\frac{\alpha +1}{2}}j^{\frac{\alpha 1}{2}}.$$
(22)
Again we find that the respective choices $`\alpha <1`$ and $`\alpha >1`$ lead to monotonically decreasing and increasing numbers of links per node as a function of node number, while setting $`\alpha =1`$ ensures the number of links per node is independent of node number. In this case, the preferential attachment of links to earlier nodes does indeed act to cancel the initial bias in link number towards later nodes. It is also apparent that when $`\alpha =1`$, the limit $`l1`$ implies all nodes possess exactly $`N`$ links as expected for a fully connected regular network. (Preferential attachment cannot distort connectivity numbers in this case as all nodes have an equal number of links.) Additionally, in the limit $`l0`$ we have $`t_{jN}=0`$ as required for an entirely disconnected network. The case $`\alpha =0`$ duplicates results found for growing networks which add a constant number of links with each new node subject to preferential attachment Albert\_02\_47 .

As previously, it is straightforward to calculate the final outbound link distribution in the case $`\alpha =1`$ using Eq. 11. This gives

$`T(k,N)`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle _0^N}𝑑j\delta (klN)`$ (23)
$`=`$ $`\delta (klN).`$
Again, we find the expected compact distribution resulting when all nodes possess the same average number of links. This raises the question however, of how it is that a probabilistic accelerating network subject to preferential attachment can end up with all nodes possessing the same average number of links? The answer lies in our use of two classes of distinguishable nodes, regulators and non-regulators, which requires that we take into account the known distribution of regulators with node number over the genome. The average link number per node at node $`n_j`$ (Eq. 22) equates to the product of the average number of link tails per regulator at node $`n_j`$, denoted $`t_r(j,N)`$, and the average number of regulators per node at node $`n_j`$, denoted $`\rho (j)`$. This latter density is $`\rho (j)=dR(j)/djlj`$ by Eq. 6, so by definition, we have

$$t_{jN}=t_r(j,N)\rho (j),$$
(24)
giving

$$t_r(j,N)=N^{\frac{\alpha +1}{2}}j^{\frac{\alpha 3}{2}}.$$
(25)
Hence, for $`\alpha <3`$, the average number of links per regulator is a decreasing function of node number $`j`$ as the growing number of links added to recent nodes is insufficient to outweigh the effects of preferential attachment which more rapidly increases the number of links attached to early nodes. In particular, for $`\alpha =1`$ with the addition of a linearly increasing number of links per node, the average number of regulatory links per regulator scales inversely with node number $`j`$. In other words, the density of regulators is very low at small node numbers $`j`$ while the very few regulatory nodes in this stretch of the genome are heavily connected due to preferential attachment so as to maintain the constant average of Eq. 22. (See Fig. 2.)

The $`t_r(j,N)`$ distribution contains information about both node connectivity and node age and so approximates genome statistics (simulated or observed) when all of this information is available. However, it is usually the case that node age information is unavailable necessitating calculation of connectivity distributions that are not conditioned on node age. This effectively requires binning together all nodes irrespective of their age to obtain a final link distribution. In the case of linearly growing number of links per node, $`\alpha =1`$, the delta function of Eq. 11 is resolved by the equality $`j=N/k`$ giving the final distribution as

$`T(k,N)`$ $`=`$ $`{\displaystyle \frac{1}{N}}\left({\displaystyle \frac{j}{k}}\right)`$ (26)
$`=`$ $`{\displaystyle \frac{1}{k^2}},`$
which, as required, is normalized as $`_1^{\mathrm{}}T(k,N)=1`$. The expected proportion of regulators $`P_t(k)`$ possessing $`k`$ links is then obtained by integrating the continuous distribution of Eq. 26 over appropriate ranges $`[1,3/2]`$ or $`[k1/2,k+1/2]`$ to obtain

$$P_t(k)=\{\begin{array}{cc}\frac{1}{3}& k=1\\ & \\ \frac{4}{4k^21}& k>1.\end{array}$$
(27)
These theoretical predictions compare well to simulations of networks of various sizes with linearly increasing numbers of probable links per node and subject to preferential attachment. Fig. 5 shows simulated outbound link distributions which are long-tailed and scale free with probabilities scaling roughly as $`P_t(k)k^2`$ for large $`k`$. The $`P_t(k)`$ distribution shows a full one third of regulators have only one link, while 60% have two or fewer links, and 71% have three or fewer links. Fig. 6 shows the long-tailed distribution $`P_t(k)`$ expected for a simulated E. coli network of $`N=2528`$ operons with preferential attachment of links. This figure shows marked similarities to the long-tailed distribution of E. coli shown in Fig. 2(c) of Ref. Cherry\_03\_40 . In particular, the expected number of regulators with $`k`$ links is $`P_t(k)R(N)`$ with the number of regulators $`R(N)`$ obtained from Eq. 6 (or from observation). For E. coli, this predicts the probable existence of about one regulator possessing link numbers in each of the respective ranges between $`[40,49]`$ links, between $`[50,64]`$ links, between $`[65,94]`$ links, between $`[95,169]`$ links, and between $`[170,700]`$ links for instance. (This approximates the connectivity of the global food sensor CRP which regulates up to 197 genes directly MandanBabu\_03\_1234 .) The average of the $`P(k)`$ distribution (as well as the $`t_r(j,N)`$ distribution) is formally undefined as long as the integration limits are taken to infinity. However, in a network of $`N`$ nodes, a regulator can practically only regulate a total of $`N`$ nodes, and this cutoff allows us estimate the average connectivity per regulator (complementing previous estimates following Eq. 7). Using the cutoff and approximating the summation via an integral, the average connectivity per regulator in a network of $`N`$ nodes is

$`k`$ $`=`$ $`{\displaystyle \underset{k=1}{\overset{N}{}}}kP_t(K)`$ (28)
$`=`$ $`{\displaystyle \frac{1}{3}}+{\displaystyle \frac{1}{2}}\mathrm{ln}\left({\displaystyle \frac{4N^21}{15}}\right),`$
(or simply $`\mathrm{ln}N`$ using the continuous distribution of Eq. 26.) The average number of links per regulator for E. coli from Eq. 28 is $`k=7.51`$ (or 7.83 using the simpler derivation), which again compares well to the observed value of 5 in E. coli Shen\_Orr\_02\_64 .

## IV Inherent prokaryote size limits

The accelerating nature of regulatory gene networks necessarily means that these networks must exhibit a transition at some critical network size either to a nonaccelerating architecture permitting continued growth or must cease growth entirely, and we now seek to predict the location of this transition point and compare it to the evolutionary record. We begin by examining an overview of the accelerating genome model. Fig. 7 shows that linear growth in link numbers per node ($`\alpha =1`$) allows a quadratic growth in the total number of links (Eq. 4) despite each of the number of regulators (Eq. 6) and the number of regulated nodes (Eq. 15) asymptoting to some fraction of $`N`$ after an initial period of quadratic growth. For large genomes, almost all new nodes will be regulators and densely connected into the existing network which will then multiply regulate almost every node.

The transition from sparse to dense connectivity occurs as an increasing proportion of operons become linked into the regulatory network leading to the emergence of a single giant component of fully connected nodes. One way to highlight this transition is by determining the proportion of transcription factor which control downstream regulators as such linkages create the single giant component. The proportion of regulators controlling regulators is

$`P_{rr}(N)`$ $`=`$ $`{\displaystyle \frac{1}{R(N)}}{\displaystyle \underset{k=1}{\overset{N}{}}}\left[1(1l)^k\right]{\displaystyle \frac{N}{k}}{\displaystyle \frac{R(N)}{N}}`$ (29)
$``$ $`lN.`$
Here, the first fraction on the RHS normalizes the proportion in terms of the number of regulators $`R(N)`$ (Eq. 6), the first term in the summation is the probability that node $`n_k`$ is a regulator, the second term is the average number of regulatory outbound links for this regulatory node $`t_r(k,N)`$ at network size $`N`$ (Eq. 25 with $`\alpha =1`$), and the third term approximates the probability that these nodes link to one of the existing regulators under random attachment. (If the very first and very last terms are dropped, the remaining summation over all nodes of the probability that $`n_k`$ is regulatory with the stated number of links equates to the total number of links in the network $`LlN^2`$. This is the more accurate version of the calculation leading to Eq. 25.) Hence, the proportion of regulators which control transcription factors scales linearly with network size and equals 15% for an $`N=2000`$ network, 29% for $`N=4000`$, 44% for $`N=6000`$, 59% for $`N=8000`$, 73% for $`N=10000`$, and 88% for $`N=12000`$ operons (after which the approximations made break down). Naturally, when most regulators themselves control other regulators, then the entire regulatory network will consist of a single giant component. These ratios compare reasonably well with those observed in E. coli where Ref. MandanBabu\_03\_1234  noted that of 121 transcription factors for which one or more regulatory genes are known, 38 factors or 31.4% regulate other transcription factors. The approximate second line of Eq. 29 with $`N=2528`$ for E. coli determines this proportion as $`P_{rr}=18.5\%`$ while the more accurate top line gives the proportion of regulators which control transcription factors as $`P_{rr}=17.7\%`$, giving a reasonable match between prediction and observation.

As the proportion of regulators of transcription factors rises, the probable length of regulatory cascades will increase. In fact, the proportion of regulators taking part in a regulatory cascade of length $`n1`$ is

$$p_n=(1P_{rr})P_{rr}^{n1}.$$
(30)
This equation can be obtained from a tree of all binary pathways which at each branching point either terminate with probability $`(1P_{rr})`$ or cascade with probability $`P_{rr}`$. As such, the probable cascade length is negligible when the proportion of regulators controlling regulators is small $`P_{rr}1`$ but can become large as $`P_{rr}`$ itself increases. As $`P_{rr}`$ is indeed large for networks of size $`N>6000`$, this again suggests that long cascades of regulatory interactions will lead to the coalescing of a single giant component in this regime. Again, the calculated lengths of regulatory cascades can be compared to those in E. coli where the number of cascades of regulated transcription factors observed in a particular set of regulatory interactions was 23 two-level cascades or 37.7%, 32 three-level cascades or 52.5%, and 6 four-level cascades or 9.8% MandanBabu\_03\_1234 . As one-level or autoregulatory interactions are not included in this observation, the predicted proportions for E. coli are $`\overline{p}_n=p_n/(1p_1)`$ with $`P_{rr}=17.7\%`$ giving 82% two-level cascades, 15% three-level cascades, 3% four-level cascades, 1% five-level cascades, and so on. It is seen that the theoretical predictions overestimate the proportion of two-level cascades and underestimate the number of three-level and higher cascades probably because of selection pressures not included in the model. Lastly, we note that the number of cycles involving closed regulatory loops of size greater than one (i.e. involving more than autoregulation) in the examined portion of the E. coli regulatory network is zero reflecting that feedback loops in these organisms are carried out at the post-transcriptional level involving metabolites such as appear in the lac operon Thieffry\_98\_43 ; Shen\_Orr\_02\_64 ; MandanBabu\_03\_1234 .

We note that our model is entirely unable to explain the high proportion of autoregulation observed in E. coli with various estimates that 28.1% Rosenfeld\_02\_785 , 50% Shen\_Orr\_02\_64  and 46.9% Thieffry\_98\_43  of regulators are autoregulatory. The predicted proportion of autoregulators is approximated by replacing the very last fraction ($`R/N`$) in Eq. 29 by the term $`1/N`$ giving the probability that a self-directed link is formed, leading to the expected autoregulatory proportion $`2/N0.08\%`$ for E. coli. This failure likely reflects the action of selection processes promoting spatial rearrangements of entire regulons on the genome and the internal shuffling of genes and promotor units. Such reorganizations of duplicated gene regions (presumably shuffling genes and promotor regions) have been common in E. coli allowing for instance, spatial regulatory motifs whereby the promotors of colocated (overlapping) and often co-functional operons transcribed in opposing directions can interfere Warren\_03\_0310029 .

The transition point from sparse to dense connectivity can be roughly located using the continuum approximation Barabasi\_99\_17 ; Barabasi\_99\_50 ; Dorogovtsev\_01\_056125 . These methods have not previously been used for this purpose (to our knowledge) and we first validate their use by deriving the known result that non-growing random graphs of $`N`$ nodes connected by an increasing number of $`L`$ undirected links undergo a phase transition from sparse to dense connectivity when $`L=N/2`$ Erdos\_60\_17 . As the number of links $`L`$ grows, the $`N`$ nodes are interlinked into firstly separate islands of size $`s_i`$ nodes for $`i=1,2,\mathrm{}`$ which eventually link up to form a giant component designated $`s_1`$ containing essentially all nodes $`s_1N`$. The largest component grows whenever a newly added link has either its head or tail in island $`s_1`$ (with probability $`s_1/N`$) and the other outside it (with probability $`(Ns_1)/N`$) leading to a size increment equal to the average size of the external islands ($`s_{j1}`$), giving

$$\frac{ds_1}{dL}=\left[2\frac{s_1}{N}\frac{(Ns_1)}{N}\right]s_{j1}.$$
(31)
Numerical or analytic integration of this equation with initial conditions $`s_1=2`$ when $`L=1`$ and assuming the average size of external smaller islands is $`s_1/2`$ shows the largest island saturating the entire network when $`L=N/2`$ as expected. (This simple approach is indicative only and is quite sensitive to for instance, the assumed average size of external islands.)

This result suggests the following transition point in directed regulatory gene networks. Each undirected (i.e. bidirectional) link in random graph theory is equivalent to two directed links allowing bidirectional traffic between any two nodes, suggesting a transition point in directed graphs at roughly $`L=N`$. This analysis suggests that the largest component is expected to saturate the entire network when link number $`LN`$ or $`N=1/l=13,677`$ (see Fig. 7). In turn, this suggests that for $`N<13,677`$ a typical network likely consists of isolated trees, while if $`N>13,677`$ the network likely consists of a single giant cluster where almost every node is connected to all others via intermediate links. When the link number is very large, $`N13,677`$, then the network becomes regularly connected Albert\_02\_47 . As prokaryote regulatory networks likely consist of functionally distinct regulated modules Thieffry\_98\_43 ; Hartwell\_99\_c4 , it is unlikely that prokaryotic gene networks can successfully operate in the fully connected regime suggesting that prokaryote genome sizes are size constrained $`N13,677`$. In fact, the previously noted absence of regulatory cycles in E. coli Thieffry\_98\_43 ; Shen\_Orr\_02\_64 ; MandanBabu\_03\_1234  likely reflects the evolutionary importance of maintaining disjoint and non-interfering regulatory units.

These results of random graph theory are suggestive only, and we now turn to consider the size of the largest connected island in prokaryote gene networks featuring directed links whose tails are preferentially attached to regulators and whose heads are randomly distributed over all existing nodes. A further difference is that prokaryote regulatory networks are themselves growing with each added node accompanied by a probabilistic number of links. In addition, we define an island to consist of all nodes which are linked regardless of the orientation of all links and so effectively treat links as being undirected. This is because a regulator can potentially perturb every node downstream to it including those nodes downstream of other regulators and so can modify the regulatory effects of other regulators—essentially, if the downstream effects of different regulators eventually intersect, we count these regulators in the same island. (Other definitions of islands could be used.)

The dominant (but not sole) mechanism by which island $`s_1`$ can grow is for the newly added node $`n_k`$ to either (a) be a regulator (with probability $`[1(1l)^k]`$) and establish an outbound regulatory link to some existing node in $`s_1`$ (with probability $`s_1/k`$) while at the same time accepting a regulatory link (with probability $`[1(1l)^k]`$) from a node in a different island $`s_{j1}`$ (with probability $`(ks_1)/k`$), or (b) accept an inbound regulatory link (with probability $`[1(1l)^k]`$) from a regulator in island $`s_1`$ (with probability $`s_1/k`$) while establishing a regulatory link (with probability $`[1(1l)^k]`$) to some node in a different island $`s_{j1}`$ (with probability $`(ks_1)/k`$). (Here, we assume that regulators are uniformly distributed over islands and the number of links within an island scales with the size of the island to crudely model preferential attachment.) The result is that island $`s_1`$ grows by the size of the second island assumed to be $`s_{j1}`$. Altogether, the rate of growth in the size of island $`s_1`$ is then

$$\frac{ds_1}{dk}=2\left[1(1l)^k\right]^2\frac{s_1[ks_1]}{k^2}s_{j1}.$$
(32)
For initial conditions, we assume that a first link appears when the genome has $`(pg_0^2)^{1/2}=177`$ nodes ($`s_1(177)=2`$). Simulations show that sufficient small islands are created to ensure $`s_{j1}`$ remains roughly constant and equal to $`s_{j1}=2.72`$, though matching the simulated and predicted curves at the 50% point requires setting $`s_{j1}=30`$. This is reasonable given the approximations made. Fig. 8 shows the size of the largest island $`s_1`$ as a proportion of all nodes. A single giant component is expected to form at a critical genome size of $`N_c=9,029`$ operons defined as the point where the simulated proportion of nodes in the giant component is 50%. (Choosing a parameter setting of 40% would also be justifiable and would lead to an exact match between predicted and observed maxima.) Unlike random graph theory, this critical point applies to all growing genomes as it is determined by the value of the link formation probability $`l`$. Genomes of smaller size than this critical value $`N<N_c`$ are expected to be sparsely connected so the network consists of multiple discrete connected islands (as in E. coli Thieffry\_98\_43 ), while genomes of larger size $`N>N_c`$ are expected to be densely connected into a single giant component where every regulator eventually perturbs the downstream effects of every other regulator.

Simulations of example genomes of various sizes spanning this critical network size confirm the adequacy of the continuum treatment. Fig. 8 shows the number of all discrete islands as well as the number of islands containing two, three and four components. In the vicinity of the critical genome size $`N_c=9,029`$, the number of discrete interconnected islands begins to decline as the growing number of links connects more and more islands into the single giant component. The size of the simulated giant component as a proportion of genome size is also shown. This figure suggests that the E. coli genome of $`N=2528`$ operons should possess a giant component containing about 5% of all nodes (about 100 nodes) which can be compared with the observation that about 70% or 300 operons of the examined regulatory and regulated operons (but not including unregulated and nonregulatory operons) could be loosely grouped into 3-6 “dense overlapping regulons” or DORS while the remaining operons appeared as disjoint systems with most containing 1-3 operons but some containing up to 25 operons Shen\_Orr\_02\_64 .

The critical network size of $`N_c=9,029`$ operons or about $`N_g=15,349`$ genes corresponds to the point where growing regulatory networks exploiting accelerating links can no longer maintain discrete functional units, islands, of interconnected nodes. Larger genomes are densely connected into a single giant component where, eventually, any regulator can perturb the downstream effects of every other node so for instance, it is unlikely that the discrete network motifs found in the E. coli regulatory network Shen\_Orr\_02\_64  can survive in this regime. This massive increase in perturbative effects immeasurably increases the difficulty of the evolutionary search process, leading to an expectation that the rate of evolutionary change will drastically slow when growing genome sizes reach criticality $`NN_c`$. From a biological point of view, it is relatively easy to understand why the critical network size $`N_c`$ acts as an upper size limit. The accelerating nature of the prokaryote regulation network means that larger networks can add new nodes only be integrating an increasing number of links to gain evolutionary benefits. Of course, the probability of finding $`lN`$ beneficial links is a rapidly decreasing function of $`N`$. It is relatively easy to find a beneficial regulator making only of order one link to existing genes (only billions of trials are needed say), but much harder when the regulator is making an average of five links with existing genes (many trillions of trails are needed). Essentially, the more links that must be beneficially integrated, the longer the evolutionary search task and the slower the rate of evolutionary change.

Many other statistical measures suggest that the regulatory mechanisms optimized to perform in a sparsely connected network will not necessarily operate in a densely connected network—evolution cannot foresee later needs. In particular, the proportion of operons $`n_j`$ which are regulated by $`k`$ inputs is, using Eq. 13, given by

$$P_o(k,N)=\frac{1}{N}\underset{j=1}{\overset{N}{}}H(k,N)=\left(\genfrac{}{}{0pt}{}{N}{k}\right)l^k(1l)^{Nk}.$$
(33)
This distribution increases with increasing network size and is shown in Fig. 9 making it clear that small networks mainly possess operons which are either entirely unregulated or regulated by only one or a few regulators. In contrast, large networks ($`N>N_c`$) have only a small proportion of operons which are unregulated while the majority of operons take between one or more regulatory inputs. It is a more difficult evolutionary task to integrate many inputs to achieve a beneficially regulated output again suggesting that prokaryote regulatory networks featuring accelerating growth in link number are size limited due to their regulatory architecture.

Another way to suggest the strict size limits imposed by the accelerating growth of regulatory links is to consider the probability that the most recently added node $`n_N`$ in a network of size $`N`$ immediately becomes regulatory. Using Eq. 5, node $`n_N`$ is a regulator with probability

$$P_r(N)=\underset{i=1}{\overset{N}{}}P(i,N)=\mathrm{\hspace{0.33em}1}(1l)^N.$$
(34)
This probability tends to unity as network size increases, and in particular, surpasses about 50% when networks consist of $`N_c`$ operons—see Fig. 10. At about this stage, large networks cannot add a new node without it having a significant probability of modifying the dynamics of existing nodes. This immeasurably increases the difficulty of the evolutionary task and again suggests a maximum size limit to prokaryote gene regulatory networks.

If the accelerating regulatory networks of prokaryotes were able to operate in the densely connected regime, the evolutionary record might be expected to show prokaryotes of arbitrarily large genome size with a transition in connectivity statistics at some critical genome size of about $`N_c9,029`$. Conversely, should these regulatory networks be unable to operate in the densely connected regime, then the evolutionary record should show a maximum size limit to prokaryote genome sizes of about $`N_c9,029`$ operons or about $`N_g=15,349`$ genes, close to the observed upper limit.

## V Conclusion

In this paper, we generalize models of accelerating networks by including probabilistic links to allow arbitrarily rapid acceleration rates leading to structural transitions in growing networks sometimes severe enough to strictly constrain network size. These structural transitions from sparse to dense connectivity are made more difficult by any additional steric or logical limitations on combinatoric control at any given promotor. Such transitions are in sharp contrast to the stationary statistics and unbounded growth potential of non-accelerating scale free and exponential networks. These probabilistic accelerating networks were applied to model prokaryote regulatory networks which exploit a quadratic growth in the number of regulators and regulatory links with genome size as established via comparative genomics programs. Our models predict a maximum genome size of $`N_c9,029`$ operons or about $`N_g=15,349`$ genes for prokaryotes, closely approximating the observed maximum. We further validated our model by making a detailed comparison of predicted and observed results for E. coli, and achieved satisfactory matches for respectively, the number of observed regulators, an average promotor binding site length of about 7, the long tailed distribution of outgoing regulatory links with an average of between 2.12 and 7.51 (compared to 5), the exponential distribution of incoming regulatory links with an average of around 1.10 (compared to 1.5), the proportion of regulators controlling regulators of around 17.7% (compared to 31.4%), and the probable length of regulatory cascades and the absence of regulatory loops. Our approach is unable to explain the high proportion of autoregulation observed in E. coli Shen\_Orr\_02\_64  and this failure likely points to selection for genome reorganizations leading to spatial arrangements of operons allowing joint regulation Warren\_03\_0310029  which is not included in this model. Further, this approach does not include selection pressures ensuring that similarly regulated islands or modules share common functionality Shen\_Orr\_02\_64 , or other regulatory mechanisms influencing both the transcription and translation of transcription factors including micro-RNAs and other chemical mechanisms and mediators (see for instance Vogel\_03\_6435 ).

However, the many successes of the accelerating network model of prokaryote regulatory networks are meaningless if similar results can be achieved via non-accelerating network models. In later work, we will show that the two simplest non-accelerating network models fail to explain either the observed quadratic growth of regulator number with genome size or the detailed statistics pertaining to the E. coli genome Gagen\_0312022 . In addition, the simplifying assumption adopted here that gene duplications ensure that operons become regulatory only on entry to the genome will be dropped in later work. This will develop a more realistic model including separate physical processes for transcription factor binding to DNA and for establishing regulatory links with regulated operons where links can form at any time.

This work has wider significance due to the still common presumption in molecular biology that “What was true for E. coli would also be true for the elephant” capturing the notion that the mechanisms operating in prokaryotes are exactly identical to those operating in complex multicellular eukaryotes. In this picture, eukaryotes are merely enlarged prokaryotes. The results of this paper indicate that this is not possible—the accelerating nature of regulatory networks necessarily implies that eukaryotes cannot be scaled up prokaryotes and that the (likely) accelerating regulatory networks of eukaryotes must be exploiting novel regulatory mechanisms. The successful modelling of these mechanisms will likely require incorporating computationally complex technologies Mattick\_01\_1611 ; Mattick\_01\_986 ; Mattick\_03\_930  into an accelerating network model, and this also will be addressed in later work.