text:
# Small world effects in evolution

## I Introduction

Darwinian evolution on asexual organisms acts by two mechanisms: mutations, that increase the genetic diversity, and selection, that fixes and reduces this diversity.

We classify mutations into point ones, corresponding to local displacements on the genotypic space (defined more accurately in the following), and other type of mutations or rearrangements, which imply larger jumps. Mutations are quite rare, so we can safely assume that for each generation at most only one mutation occurs. From a biological point of view, this discussion applies to unicellular asexual organisms, for which there is no distinction between somatic and germ cells. Moreover, we do not consider the possibility of recombination (exchange of genetic material). For non-recombinant (asexual) organisms, the combined effects of reproduction and mutations correspond to a random walk on the genotypic space. Even sexual (recombinant) populations do transmit asexually some part of their genotype, such as mithocondrial or y-chromosome DNA, to which the following analysis applies. Furthermore, we assume a constant environment.

We schematize the effects of selection by the selective reduction of the survival probability. Selection results from several constraints and can affect the frequency of certain genotypes, or even eliminate them from the population. We classify the components of selection into static and dynamic ones. In the first class we put the constraints that are independent on the population distribution, such as the functionality of a certain protein or the tuning of a metabolic path. The static part of the selection is equivalent to the concept of “fitness landscape” . The dynamic part of the selection originates from to the competition among individuals. We assume in this schematization that the competition only arises because the population shares some external resource which is evenly distributed, i.e. we disregard the effects of the inter-individual competition . In this limit the effects of competition do not depend on the distribution of genotypes, and simply limit the total size of the population.

The most important effects of the selection can be roughly schematized by assuming that certain genotypes are forbidden, so that they are eliminated from the accessible space. Let us consider for the moment that evolution takes place in a flat fitness landscape. In this frame, there are no interactions among individuals: the evolution is given by the simple superposition of all possible lineages. The probability distribution of the population in the genotypic space at a given time can be obtained by summing over all possible individual “histories” in a way similar to the path integral formulation of statistical mechanics . In this view, the phylogenetic lineage of an individual is given by the path connecting the genotypes of ancestors.

This assumption has to fail somewhere, since otherwise evolution would be equivalent to a diffusion process, without anything favoring the formation of species. It often assumed that speciation events are rare and occur in a very short time , either due to some change in the fitness landscape (caused by an external catastrophes like the fall of a large meteorite, or by an internal rearrangement induced by coevolution) or because small isolated populations, escaping from competition, are free to explore their genotypic space (genetic drift) and find a path to an higher fitness peak . Therefore, it is usually assumed that one can reconstruct the diverging time of two species from their genotypic distance from a single speciation event .

It is however known from recent investigations about the small-world phenomenon , that diffusion is seriously affected even by a small fraction of long-range jumps . This fact could have dramatic consequences in our understanding of the large-scale evolution mechanism. If the short-range mutations are dominant, the evolution is equivalent to a diffusion process in the genotypic space. Assuming that the fitness landscape is formed by mountains separated by valleys, and that the crests of mountains are almost flat, then the large-scale evolution is dominated by the times needed to cross a valley by chance, while the short-scale is dominated by the neutral exploration of crests .

Vice versa, if the long-range mutations are important (eventually amplified by the small-world mechanism), the time needed to connect two genotypes does not depend on their distance nor by the shape of the fitness landscape and the fitness maxima are quickly populated. Moreover, in this scenario, the speciation phenomenon should not be ascribed to the “discovery” of a preexisting niche, but rather to the formation of that niche in a given ecosystem due to internal interactions (coevolution) or external physical changes (catastrophes). After the formation, the niche is quickly populated because of the long-range mutations. In this framework, allopatric speciation  looses its fundamental importance, and sympatric speciation due to coevolution  becomes a plausible alternative.

In order to evaluate quantitatively the relevance of these speculations, we introduce an individual-based model of an evolving population. We assume that the genetic information (the genotype) of an individual is represented by a binary string $`𝒙=(x_1,x_2,\mathrm{},x_L)`$ of $`L`$ symbols $`x_i=0,1`$ (multilocus model with two alleles). In this way we are modeling haploid organisms, i.e. bacteria or viruses or, more appropriately, more archaic, pre-biotic entities. The choice of a binary code is not fundamental but certainly makes things easier. It can be justified by thinking at a purine-pyrimidine coding, or to “good” and “bad” alleles for genes (units of genetic information). In this second version, $`0`$ represents a good gene and $`1`$ a bad one. This is a sort of “minimal model” which is often used to model the evolution of genetic populations . For bacteria, one can give indications that major and minor codons work in such a way . The genotypic space is thus a Boolean hypercube of $`L`$ dimensions, and each individuals sits on a corner of this cube, according to its genotype. The $`2^L`$ corners of this hypercube represent all possible genotypes, which are at a maximum distance equal to $`L`$.

In the next section, we formalize the model and introduce the mathematical representation of mutation mechanisms. Then, in Section III, we consider the consequences of short and long-range mutations for a flat fitness landscape. In the long-range case all genotypes are connected, regardless of their mutual distance: this case can be considered the equivalent of a mean-field approximation. We derive analytically the rate of spreading $`v`$ and the characteristic spreading time $`\tau `$ of a genetically homogeneous inoculum in the short-range ($`v_s`$, $`\tau _s`$) and long-range ($`v_{\mathrm{}}`$, $`\tau _{\mathrm{}}`$) cases. In the first case, $`v_s`$ is independent of the genotype length $`L`$, and $`\tau _s`$ grows linearly with $`L`$; the opposite happens in the long-range case. Clearly, this different behavior has dramatic evolutionary consequences as $`L`$ becomes large.

In general, however, only a small set of all possible long-range mutations is observed in real organisms. We thus compute numerically the rate of spreading $`v`$ in a mixed short-range and sparse long-range case. A small-world effect can be observed: in the limit of large genotype lengths, a vanishing fraction of long-range mutations cooperates with the short-range ones to give essentially the mean-field results.

In Section IV, we show that the effects of mutations and selection can be separated in the limit of a very smooth fitness landscape and mean-field long-range mutations. In this approximation we obtain analytically that the asymptotic probability distribution is a Boltzmann equilibrium one, in which the fitness plays the role of of an energy and mutations correspond to temperature. We compute numerically the asymptotic probability distribution for several fitness landscapes, finding that the equilibrium hypothesis is verified. As a nontrivial example of small-world effects in evolution, we checked numerically that this scenario still holds for a sparse long-range mutation matrix. Conclusions and perspectives are drawn in the last section.

## II The model

In order to be specific, let us assume that the genotype of an individual is represented by a Boolean string of length $`L`$. In this way, the genotypic space is a Boolean hypercube with $`2^L`$ nodes.

The genotype of an individual is represented by a string $`𝒙=(x_1,x_2,\mathrm{},x_L)`$ of $`L`$ Boolean symbols $`x_i=0,1`$. Each position corresponds to a locus whose gene has two allelic forms. In this way we are modeling haploid (only one copy of each gene) organisms, i.e. bacteria or viruses or, more appropriately, more archaic, pre-biotic entities. The genotype can be also viewed as a spin configuration. In this case we use the symbol $`\sigma _i=2x_i1`$.

We consider two kinds of mutations: point mutations, that interchange a $`0`$ with a $`1`$, and all other mutations that do not alter the length of the genotype (transposition and inversions). We define the distance $`d(𝒙,𝒚)`$ between two genotypes $`𝒙`$ and $`𝒚`$ as the minimal number of point mutations needed to pass from $`𝒙`$ to $`𝒚`$ (Hamming distance)

$`d(𝒙,𝒚)={\displaystyle \underset{i=1}{\overset{L}{}}}(x_iy_i)^2.`$

All possible genotypes of length $`L`$ are distributed on the $`2^L`$ vertices of an hypercube. A point mutation corresponds to a unit displacement on that hypercube (short-range jump).

The occurrence of point mutations in real organisms depends on the identity of the symbol and on its position on the genotype; in the present approximation, however, we assume that all point mutations are equally likely. Moreover, since the probability of observing a mutation is quite small, we impose that at most one mutation is possible in one generation. The probability to observe a point mutation from genotype $`𝒚`$ to genotype $`𝒙`$ is given by the short-range mutation matrix $`M_s(𝒙,𝒚)`$. By denoting with $`\mu _s`$ the probability of a point mutation per generation, we have

$$M_s(𝒙,𝒚)=\{\begin{array}{cc}1\mu _s\hfill & \text{if }𝒙=𝒚\text{,}\hfill \\ \frac{\mu _s}{L}\hfill & \text{if }d(𝒙,𝒚)=1\text{,}\hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}$$
(1)
Other types of mutations correspond to long-range jumps in the genotypic space. The simplest approximation consists in assuming all mutations equiprobable. Let us denote with $`\mu _{\mathrm{}}`$ the probability per generation of this kind of mutations. The long-range mutation matrix, $`M_{\mathrm{}}(𝒙,𝒚)`$, is defined as

$$M_{\mathrm{}}(𝒙,𝒚)=\{\begin{array}{cc}1\mu _{\mathrm{}}\hfill & \text{if }𝒙=𝒚\text{,}\hfill \\ \frac{\mu _{\mathrm{}}}{2^L1}\hfill & \text{otherwise.}\hfill \end{array}$$
(2)
In the real world, only some kinds of mutations are actually observed. We model this fact by replacing $`𝑴_{\mathrm{}}`$ with a sparse matrix $`\widehat{𝑴}_{\mathrm{}}`$. We introduce the sparseness index $`s`$, which is the average number of nonzero off-diagonal elements of $`\widehat{𝑴}_{\mathrm{}}`$. The sum of these off-diagonal elements still gives $`\mu _{\mathrm{}}`$. In this case $`\widehat{𝑴}_{\mathrm{}}`$ is a quenched sparse matrix, and $`𝑴_{\mathrm{}}`$ can be considered the average of the annealed version.

After considering both types of mutations, the overall mutation matrix is $`𝑴=𝑴_{\mathrm{}}𝑴_s`$ or $`𝑴=\widehat{𝑴}_{\mathrm{}}𝑴_s`$ for the quenched version.

We model our population at the level of the probability distribution of genotypes $`𝒑𝒑(t)`$, thus disregarding spatial effects. The evolution equation for $`𝒑`$ is

$$p^{}(𝒙)=\frac{A(𝒙)_𝒚M(𝒙,𝒚)p(𝒚)}{\overline{A}},$$
(3)
where the selection function $`A(𝒙)`$ corresponds to the average reproduction rate of individuals with genotype $`𝒙`$, and $`\overline{A}=_𝒙A(𝒙)p(𝒙)`$ is the average reproduction rate of the population . We write $`A(𝒙)`$ in an exponential form $`A(𝒙)=\mathrm{exp}(V(𝒙))`$, and we denote $`V(𝒙)`$ with the term *fitness landscape*.

The selection does not act directly on the genotype, but rather on the phenotype (how an individual appears to others). The phenotype of a given genotype can be interpreted as an array of morphological characteristics. We consider the simplest case, in which the phenotype if univocally determined by the genotype $`𝒙`$, which is not the general case, since polymorphism or age dependence are usually present. The general mapping between genotype and phenotype is largely unknown and is expected to be quite complex. The effects of some genes are additive (non-epistatic), while others can interact in a simple (control genes) or complex (morphologic genes) way.

A possible way of approximating these effects is to use the following form for the fitness $`V(𝒙)`$:

$$V(𝒙)=\frac{}{L}\underset{i=1}{\overset{L}{}}\sigma _i+\frac{𝒥}{L1}\underset{i=1}{\overset{L1}{}}\sigma _i\sigma _{i+1}+𝒦\eta (𝒙),$$
(4)
where $`\eta (𝒙)`$ is a random function of $`𝒙`$, uniformly distributed between $`1`$ and $`1`$ ($`\eta (𝒙)\eta (𝒚)=\delta _{𝒙𝒚}`$). The “field” term $``$ represents a non-epistatic contribution to the fitness, in which all genes have equal weight. The “ferromagnetic” term $`𝒥`$ represents simple interactions between pairs of genes (even though in general those are not symmetric) and corresponds to a weakly rough landscape. Finally, $`𝒦`$ modulates a widely rough landscape and can be thought as an approximation of the effects of complex (spin-glass like) interactions among genes.

## III Spreading and small-world effects on a flat fitness landscape

In this section we study the case of a flat fitness landscape (no selection), i.e. $`=𝒥=𝒦=0`$. In this way we are modeling the evolution on the crest of a mountain, assuming that all deleterious mutations are immediately lethal (the part of genotype that can originate this kind of mutations is not considered), and that we can neglect the small variation of fitness along the crest. This landscape is the one usually considered in the theory of neutral evolution . Let us assume that this crest is colonized by a founder deme genotypically homogeneous (a delta peak in the genotypic distribution). We want to obtain the average time needed to populate the crest according with the different mutation schemes. Since the fitness landscape is flat, $`\overline{A}`$ is a constant.

We are interested in the behavior of $`\rho (t)`$, which is the average distance of the population from a given starting genotype $`𝒙_0`$, i.e.

$$\rho (t)=\underset{𝒙}{}p(𝒙,t)d(𝒙,𝒙_0).$$
(5)
We introduce the spreading velocity $`v`$ as

$`v={\displaystyle \frac{\rho }{t}}|_{t=0}.`$

In Appendix A we obtain the spectral properties of the mutational matrices $`𝑴_s`$ and $`𝑴_{\mathrm{}}`$. The corresponding spreading velocities in the limit of small time interval $`t`$ compared to the characteristic spreading time $`\tau _s=L/2\mu _s`$ (short-range) and $`\tau _{\mathrm{}}=1/\mu _{\mathrm{}}`$ (long-range) are $`v_s=\mu _s`$, Eq. (12), and $`v_{\mathrm{}}=L\mu _{\mathrm{}}/2`$, Eq. (13). If short and long-range mutations coexists, one has $`v=v_s+v_{\mathrm{}}`$, Eq. (14).

One can approximate the behavior of real biological systems by considering a mixture of short-range mutations, which occur with a relatively high frequency, and sparse long-range mutations, with sparseness index $`s`$.

We investigated the sparse case by numerical simulations, for some genotype lengths $`L`$. As shown in Figure 1, as soon as the sparseness index $`s>0`$, the numerical value of $`v`$ becomes very close to the mean-field one, Eq. (14). Notice that the average distance from inoculum, $`\rho (t)`$, is rather insensitive to the distribution of $`𝒑`$. The actual distribution can be quite different from the one obtained with the mean-field matrix $`𝑴_{\mathrm{}}`$.

This transition may be interpreted as an indication of a small-world effect, and that there exists a first-order transition at $`s=0`$ . However, the fluctuations of the velocity, $`\text{Var}(v)`$, appears to diverge at $`s=0`$ with an exponent $`1`$, as shown in Figure 2.

These results suggest the presence of a small-world phenomenon in evolution: the rare and sparse long-range mutations cooperate synergetically with the short-range ones to give essentially a mean-field effect. We checked that the previous results hold also for a non-flat (but smooth) fitness landscape ($`\mu _s=10^5`$, $`\mu _{\mathrm{}}=10^4`$ and $`=10^4`$, $`𝒥=0`$, $`𝒦=0`$). Again, as soon as $`s>0`$, the spreading velocity increases from the short-range values to the mean-field ones.

## IV Equilibrium properties and small-world effects on a smooth fitness landscapes

In the following we shall study the effects of the cooperation between short-range and sparse long-range mutations on the equilibrium properties of the model.

Let us study the model in the presence of a smooth static fitness landscape. In this case the fitness $`A`$ does not depend explicitly on the distribution $`𝒑`$, and Eq. (3) can be linearized by using unnormalized variables $`z(𝒙,t)`$ satisfying

$$z(𝒙,t+1)=\underset{𝒚}{}A(𝒚)M(𝒙,𝒚)z(𝒚,t),$$
(6)
with the correspondence

$`p(𝒙,t)={\displaystyle \frac{z(𝒙,t)}{_𝒚z(𝒚,t)}}.`$

In vectorial terms, Eq. (6) can be written as

$$𝒛(t+1)=𝑴𝑨𝒛(t),$$
(7)
where $`𝑴_{𝒙𝒚}=M(𝒙,𝒚)`$ and $`𝑨_{𝒙𝒚}=A(𝒙)\delta _{𝒙𝒚}`$.

When one takes into consideration only point mutations ($`𝑴𝑴_s`$), Eq. (6) can be read as the transfer matrix of a two-dimensional Ising model , for which the genotypic element $`\sigma _i(t)`$ corresponds to the spin in row $`t`$ and column $`i`$, and $`z(𝝈,t)`$ is the restricted partition function of row $`t`$. The effective Hamiltonian $`𝒱`$ (up to constant terms) of a possible genealogical history $`\{𝒙(t)\}`$ or $`\{𝝈(t)\}`$ from time $`1tT`$ is

$$𝒱=\underset{t=1}{\overset{T1}{}}\left(\gamma \underset{i=1}{\overset{L}{}}\sigma _i(t)\sigma _i(t+1)+V\left(𝒙(t)\right)\right),$$
(8)
where $`\gamma =\mathrm{ln}(\mu _s/(1\mu _s))`$.

This peculiar two-dimensional Ising model has a long-range coupling along the row (depending on the choice of the fitness function) and a ferromagnetic coupling along the time direction (for small short range mutation probability). In order to obtain the statistical properties of the system one has to sum over all possible configurations (stories), eventually selecting the right boundary conditions at time $`t=1`$. The bulk properties of Eq. (8) cannot be reduced in general to the equilibrium distribution of a one dimensional system, since the transition probabilities among rows do not obey detailed balance. Moreover, the temperature-dependent Hamiltonian (8) does not allow an easy identification between energy and selection, and temperature and mutation, what is naively expected by the biological analogy with an adaptive walk.

### A Long-range mutations

Let us first consider the long-range mutation case. Eq. (7), reformulated according to Eq. (3), corresponds to

$`𝒛(t+\tau )=(𝑨𝑴_{\mathrm{}})^\tau 𝒛(t).`$

Since it is easier to consider the effects of $`𝑨`$ and $`𝑴_{\mathrm{}}`$ separately, let us study in which limit they commute. The norm of the commutator on the asymptotic probability distribution $`𝒑`$ is

$`[𝑨𝑴_{\mathrm{}}]={\displaystyle \underset{𝒙𝒚}{}}|[𝑨𝑴_{\mathrm{}}]_{𝒙𝒚}p(𝒚)|,`$

and it is bounded by $`\mu _{\mathrm{}}c`$, where $`c=\mathrm{max}_{𝒙𝒚}|A_{𝒙𝒙}A_{𝒚𝒚}|`$. In the limit $`\mu _{\mathrm{}}c0`$ (i.e. a very smooth landscape), to first order in $`c`$ we have

$`(𝑨𝑴_{\mathrm{}})^\tau =𝑨^\tau 𝑴_{\mathrm{}}^\tau +O(\tau ^2\mu _{\mathrm{}}c)𝑨^{\tau 1}𝑴_{\mathrm{}}^{\tau 1},`$

which is the analogous of the Trotter product formula.

When $`\tau `$ is order $`1/\mu _{\mathrm{}}`$, $`𝑴_{\mathrm{}}^\tau `$ is a constant matrix with elements equal to $`1/2^L`$, and thus $`𝑴_{\mathrm{}}𝒑`$ is a constant probability distribution. If $`\mu _{\mathrm{}}`$ is large enough, $`(𝑨𝑴_{\mathrm{}})^\tau =𝑨^\tau 𝑴_{\mathrm{}}^\tau `$.

The asymptotic probability distribution $`\stackrel{~}{𝒑}`$ is thus proportional to the diagonal of $`𝑨^{1/\mu _{\mathrm{}}}`$:

$$\stackrel{~}{p}(𝒙)=C\mathrm{exp}\left(\frac{V(𝒙)}{\mu _{\mathrm{}}}\right)$$
(9)
i.e. a Boltzmann distribution with Hamiltonian $`V(𝒙)`$ and temperature $`\mu _{\mathrm{}}`$. This corresponds to the naive analogy between evolution and equilibrium statistical mechanics. In other words, the genotypic distribution is equally populated if the phenotype is the same, regardless of the genotypic distance since we used long-range mutations.

We have checked numerically this hypothesis, by iterating Eq. (3) for a time $`T`$ large enough to be sure of having reached the asymptotic state. We plotted the logarithm of the probability distribution $`\stackrel{~}{p}(𝒙)`$ versus the value of the Hamiltonian $`V(𝒙)`$. We computed the slope $`1/\mu _e`$ of the linear regression. The quantity $`\mu _e`$ is the effective “temperature” of the probability distribution according to the equilibrium hypothesis.

The results for the mean-field mutation matrix $`𝑴_{\mathrm{}}`$ are shown in Figure 3. We see that the equilibrium hypothesis is well verified in the limit $`\mu _{\mathrm{}}c`$; and that convergence is faster for a rough landscape. In all cases the effective temperature $`\mu _e`$ is close to the expected value, i.e. to the mutation rate $`\mu _{\mathrm{}}`$.

### B Short-range mutations

The above results only hold qualitatively for pure short-range mutations as shown in Figure 4. A small dispersion of points in figure implies that the genotypes can be divided into evenly-populated groups sharing approximately the same fitness. This is always the case for the additive fitness landscape ($``$ contribution), since in this case the position of symbols in the genotype is ininfluent, and the short-range mutations are able to homogenize the distribution inside a group. This is only approximately valid for the $`𝒥`$ contributions, since in this case the fitness depends on pairs of symbols, while mutations only act on single symbols. However, as shown in Figure 4-d, the homogeneous state is reached for a sufficiently high mutation probability. Finally, the homogeneous state is never reached for the very rough landscape case ($``$ contribution), even though the linear relation is satisfied in average.

In order to obtain a quantitatively correct prediction, one has to consider that the resulting slope $`\mu _e`$ is related to the second largest eigenvalue $`\lambda _1`$ of the mutation matrix by $`\mu =1\lambda _1`$. When the fitness only depends on the “external field” term $``$, in the asymptotic state the short-range mutations connects group of equal fitness. This fact is reflected by the vanishing dispersion of points in Figure 4-a and 4-b.

In all cases, we expect that in the limit of a very smooth landscape, $`\mu _e`$ tends towards the limiting value $`2\mu _s/L`$ of Eq. (11). This limit is reached faster when only the $``$ term is present. This implies that for large genotype, the effective temperature due to short-range mutations is vanishing.

In the opposite case, when the selection is strong, the application of the matrix $`𝑨`$ “rotates” the distribution $`𝒑`$ in a way which is practically random with respect to the Fourier eigenvectors of $`𝑴_s`$. Thus, the effective second eigenvalue of the mutation matrix is given by $`1\mu _s`$, obtained averaging over all the eigenvalues of Eq (11). Consequently, we obtain $`\mu _e\mu _s`$. This limit is reached faster in the case of a strongly disordered fitness landscape, i.e. in the limit of large $`𝒦`$. When only the $`𝒥`$ term is present, one observes an intermediate case.

### C Small-world effects

Let us now consider the case of a sparse long-range mutation matrix, coupled to a stronger short-range matrix, i.e. $`\mu _s\mu _{\mathrm{}}`$. We performed simulations with $`\mu _s=0.1`$ and $`\mu _{\mathrm{}}=0.01`$, $`L=6`$, 8 and 10 and varying $`s`$. The results are summarized in Table I.

The effects of the two kinds of mutations is additive, implying that, for small $`\mu _{\mathrm{}}`$, the distributions are visually similar to the short-range case, Figure 4. However, as soon as $`s`$ is greater than a threshold $`s_c`$, the effective temperature $`\mu _e`$ increases by the expected long-range contribution, $`\mu _{\mathrm{}}`$. This increment is very relevant when the effect of short-range mutations are vanishing, i.e. for smooth landscapes ($``$ and $`𝒥`$ contributions) and large genotypes. The effect is less relevant in the disordered case, ($`𝒦`$ contributions), since in this case the contribution to the effective temperature by the small-range mutations do not decrease with $`L`$.

By increasing the weigth of the sparse long-range mutations, also the probability distribution becomes similar to the long-range case, Figure 3. We performed simulations with $`\mu _s=0.1`$ and $`\mu _{\mathrm{}}=0.1`$, $`L=6,8,10`$ varying $`s`$. The results are shown in Table I. Notice that the value of $`s_c`$ (estimated visually from the plot of data) is rather insensitive to $`L`$ and the parameters of the fitness $`V(𝒙)`$.

We quantified the distance between the resulted distribution and the Boltzmann one Eq. (9) through the computation of the average square difference from linear regression, $`\chi ^2`$. In Figure 5 we show that $`\chi ^2`$ becomes very small as soon as $`s>s_c`$, which seems to vanish with $`L\mathrm{}`$.

This implies that, for large genotypes, a vanishing fraction of long-range mutations, coupled to small-range ones, are sufficient to establish the statistical mechanics analogy of selection and mutations.

## V Conclusions

We studied some simple models of asexual populations evolving on a smooth fitness landscape, in the presence of point mutations (small-range jumps in genotypic space) and other genetic rearrangements (long-range jumps).

We have computed analitically the spreading velocity of an initially homogeneous inoculum on a flat fitness landscape, for the short-range and the long-range mean-field (all mutations equiprobable) cases. Since in the real situation only a small set of the all possible mutations can occur, we have also considered the quenched version of the long-range mutation matrix. In this case we have shown that a small-world effect is present, since even a small fraction of quenched long-range jumps makes the results indistinguishable from those obtained by assuming all mutations equiprobable. These results still hold for a smooth fitness landscape.

We have investigated this issue further, studying the equilibrium properties of the system in the presence of a smooth fitness landscape. In this framework, it has been possibile to show that the equilibrium distribution is a Boltzmann one, in which the fitness plays the role of an energy, and mutations that of a temperature. We have checked numerically this result for different fitness landscapes, and a mean-field long-range mutation mechanism. As in the previous case, a small-world phenomenon appears, since similar results can be obtained using a combination of sparse long-range and short-range mutations.

## Acknowledgements

We wish to acknowledge our participation to the DOCS  group. The numerical simulations have been performed using the computational facilities of INFM PAIS 1999-G-IS-Firenze.

## Appendix: spectral properties of the mean-field mutation matrices

Both $`𝑴_s`$ and $`𝑴_{\mathrm{}}`$, Eqs. (1-2), are Markov matrices. Moreover, they are circular matrices, since the value of a given element does not depend on its absolute position but only on the distance from the diagonal. This means that their spectrum is real, and that the largest eigenvalue is $`\lambda _\mathrm{𝟎}=1`$. Since the matrices are irreducible, the corresponding eigenvector $`𝝃_\mathrm{𝟎}`$ is non-degenerate, and corresponds to the flat distribution $`\xi _\mathrm{𝟎}(𝒙)=1/2^L`$. In analogy with circular matrices in the usual space, one can obtain their complete spectrum using the analogous of Fourier transform in a Boolean hyper-cubic space. Let us define the “Boolean scalar product” $``$:

$`𝒙𝒚={\displaystyle \underset{i=1}{\overset{L}{}}}x_iy_i,`$

where the symbol $``$ represent the sum modulus two (XOR) and the multiplication can be substituted by an AND (which has the same effect on Boolean quantities). This scalar product is obviously distributive with respect to the XOR:

$`(𝒙𝒚)𝒛=(𝒙𝒛)(𝒚𝒛).`$

Note that the operation $`𝒙𝒚`$ is performed bitwise between the two genotypes: $`(𝒙𝒚)_i=x_iy_i`$.

Given a function $`f(𝒙)`$ of a Boolean quantity $`𝒙\{0,1\}^L`$, its “Boolean Fourier transform” is $`\stackrel{~}{f}(𝒌)`$ ($`𝒌\{0,1\}^L`$)

$`\stackrel{~}{f}(𝒌)={\displaystyle \frac{1}{2^L}}{\displaystyle \underset{𝒙}{}}(1)^{𝒙𝒌}f(𝒙).`$

The anti-transformation operation is determined by the definition of the Kronecker delta

$`\delta _{𝒌\mathrm{𝟎}}={\displaystyle \frac{1}{2^L}}{\displaystyle \underset{𝒙}{}}(1)^{𝒙𝒌},`$

and is given by

$`f(𝒙)={\displaystyle \underset{𝒌}{}}(1)^{𝒙𝒌}\stackrel{~}{f}(𝒌).`$

One can easily verify that the Fourier vectors $`\xi _k(𝒙)=(1)^{𝒙𝒌}`$ are eigenvectors of both $`𝑴_{\mathrm{}}`$ and $`𝑴_s`$, with eigenvalues

$$\begin{array}{cc}\hfill \lambda _\mathrm{𝟎}& =1,\hfill \\ \hfill \lambda _𝒌& =1\mu _{\mathrm{}}\frac{\mu _{\mathrm{}}}{2^L1},\hfill \end{array}$$
(10)
for the long-range case, and

$$\begin{array}{cc}\hfill \lambda _\mathrm{𝟎}& =1,\hfill \\ \hfill \lambda _𝒌& =1\frac{2\mu _sd(𝒌,\mathrm{𝟎})}{L},\hfill \end{array}$$
(11)
for the short-range case, where $`d(𝒙,𝒚)`$ is the Hamming distance between genotypes $`𝒙`$ and $`𝒚`$.

The computation of $`\rho (t)`$, Eq. (5), is easily performed in Fourier space, using the analogous of Parsifal theorem:

$`\begin{array}{cc}\hfill {\displaystyle \underset{𝒙}{}}f(𝒙)g(𝒙)& ={\displaystyle \underset{𝒙}{}}{\displaystyle \underset{𝒌}{}}{\displaystyle \underset{𝒌^{}}{}}\stackrel{~}{f}(𝒌)\stackrel{~}{g}(𝒌^{})(1)^{𝒙(𝒌𝒌^{})}\hfill \\ & =2^L{\displaystyle \underset{𝒌}{}}{\displaystyle \underset{𝒌^{}}{}}\stackrel{~}{f}(𝒌)\stackrel{~}{g}(𝒌^{})\delta _{𝒌𝒌^{}}\hfill \\ & =2^L{\displaystyle \underset{𝒌}{}}\stackrel{~}{f}(𝒌)\stackrel{~}{g}(𝒌),\hfill \end{array}`$

where we have used the property $`𝒌𝒌^{}=\mathrm{𝟎}`$ if and only if $`k_i=k_i^{}`$ for each component $`i`$.

Let us denote by $`𝒆_n`$ the unit vector along direction $`n`$, i.e. $`(𝒆_n)_i=\delta _{ni}`$.

The Fourier transform of the distance $`d(𝒙,𝒚)`$ is obtained considering that $`(1)^{(𝒙𝒚)𝒆_n}`$ gives $`1`$ if $`x_n=y_n`$ and $`1`$ otherwise, thus

$`d(𝒙,𝒚)={\displaystyle \frac{1}{2}}\left(L{\displaystyle \underset{n=0}{\overset{L1}{}}}(1)^{(𝒙𝒚)𝒆_n}\right).`$

We obtain

$`\begin{array}{cc}\hfill \stackrel{~}{d}_{𝒙_0}(𝒌)& ={\displaystyle \frac{1}{2^L}}{\displaystyle \underset{𝒙}{}}d(𝒙,𝒙_0)(1)^{𝒙𝒌}\hfill \\ & ={\displaystyle \frac{1}{2^L}}{\displaystyle \frac{1}{2}}{\displaystyle \underset{𝒙}{}}\left(L{\displaystyle \underset{n=0}{\overset{L1}{}}}(1)^{(𝒙𝒙_0)𝒆_n}\right)(1)^{𝒙𝒌}\hfill \\ & ={\displaystyle \frac{L}{2}}\delta _{𝒌\mathrm{𝟎}}{\displaystyle \frac{1}{2}}\left({\displaystyle \underset{n=0}{\overset{L1}{}}}(1)^{𝒙_0𝒆_n}\right)\left({\displaystyle \underset{𝒙}{}}(1)^{𝒙(𝒌𝒆_n)}\right)\hfill \\ & ={\displaystyle \frac{L}{2}}\delta _{𝒌\mathrm{𝟎}}{\displaystyle \frac{1}{2}}{\displaystyle \underset{n=0}{\overset{L1}{}}}(1)^{𝒙_0𝒆_n}\delta _{𝒌𝒆_n};\hfill \end{array}`$

i.e.

$`\stackrel{~}{d}_{𝒙_0}(𝒌)=\{\begin{array}{cc}L/2\hfill & \text{if }𝒌=\mathrm{𝟎},\hfill \\ 1/2\hfill & \text{if }𝒌=𝒆_n\text{ and }(𝒙_0)_n=0,\hfill \\ 1/2\hfill & \text{if }𝒌=𝒆_n\text{ and }(𝒙_0)_n=1,\hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}`$

The probability distribution $`p(𝒙,t)`$ can be expanded on the eigenvector basis $`\xi _k(𝒙)`$ of $`𝑴`$:

$`p(𝒙,0)={\displaystyle \underset{𝒌}{}}a_𝒌\xi _𝒌(𝒙)={\displaystyle \underset{𝒌}{}}a_𝒌(1)^{𝒙𝒌},`$

and

$`p(𝒙,t)=𝑴^tp(𝒙,0)={\displaystyle \underset{𝒌}{}}a_𝒌\lambda _𝒌^t\xi _𝒌(𝒙),`$

i.e.

$`\stackrel{~}{p}(𝒌,t)=a_𝒌\lambda _𝒌^t.`$

Thus

$`\begin{array}{cc}\hfill \rho (t)& ={\displaystyle \underset{𝒌}{}}\stackrel{~}{d}_{𝒙_0}(𝒌)\stackrel{~}{p}(𝒌,t)\hfill \\ & ={\displaystyle \underset{𝒌}{}}\stackrel{~}{d}_{𝒙_0}(𝒌)a_𝒌\lambda _𝒌^t\hfill \\ & ={\displaystyle \frac{L}{2}}\lambda _\mathrm{𝟎}^t+{\displaystyle \underset{n=0}{\overset{L1}{}}}{\displaystyle \frac{1}{2}}a_{𝒆_n}\lambda _{𝒆_n}^t.\hfill \end{array}`$

If at $`t=0`$ the distribution is concentrated at $`𝒙_0=\mathrm{𝟎}`$ ($`p(𝒙,0)=\delta _{𝒙\mathrm{𝟎}}`$) then $`a_𝒌=1`$ for all $`𝒌`$. In both the short and long-range cases, $`\lambda _{𝒆_n}`$ does not depend on $`n`$ (see Eqs. (10)-(11)), and thus

$`\rho (t)={\displaystyle \frac{L}{2}}\left(\lambda _\mathrm{𝟎}^t\lambda _\mathrm{𝟏}^t\right).`$

For the short-range case we have

$`\begin{array}{cc}\hfill \rho _s(t)& ={\displaystyle \frac{L}{2}}\left(1\left(1{\displaystyle \frac{2\mu _s}{L}}\right)^t\right)\hfill \\ & ={\displaystyle \frac{L}{2}}\left(1\mathrm{exp}\left(t\mathrm{ln}\left(1{\displaystyle \frac{2\mu _s}{L}}\right)\right)\right)\hfill \\ & {\displaystyle \frac{L}{2}}\left(1\mathrm{exp}\left({\displaystyle \frac{t}{\tau _s}}\right)\right),\hfill \end{array}`$

and the characteristic spreading time is $`\tau _s=L/2\mu _s`$. For $`t`$ small compared with $`\tau _s`$ ($`L\mathrm{}`$) we have

$$\rho _s(t)\mu _stv_st.$$
(12)
For long-range mutations we have

$`\begin{array}{cc}\hfill \rho _{\mathrm{}}(t)& ={\displaystyle \frac{L}{2}}\left(1\left(1\mu _{\mathrm{}}{\displaystyle \frac{\mu _{\mathrm{}}}{𝒆_nL1}}\right)^t\right)\hfill \\ & {\displaystyle \frac{L}{2}}\left(1\mathrm{exp}\left({\displaystyle \frac{t}{\tau _{\mathrm{}}}}\right)\right)\hfill \end{array}`$

with $`\tau _{\mathrm{}}=1/\mu _{\mathrm{}}`$. For $`t<<\tau _{\mathrm{}}`$ ($`\mu _{\mathrm{}}0`$) we have

$$\rho _{\mathrm{}}(t)\frac{L\mu _{\mathrm{}}}{2}tv_{\mathrm{}}t.$$
(13)
The behavior of $`\rho (t)`$ for short times, vanishing mutation probabilities and large genotypes, Eqs. (12)-(13), is rather trivial. In these approximations one can neglect back mutations, and obtain the same results on a Cayley tree. However, the full analysis gives the exact behavior of $`\rho (t)`$ for all times.

In the mixed case one has $`𝑴=𝑴_s𝑴_{\mathrm{}}`$. Since $`𝑴_s`$ and $`𝑴_{\mathrm{}}`$ share the same eigenvectors, $`\lambda =\lambda _s\lambda _{\mathrm{}}`$ and, in the previous approximations,

$$v=v_s+v_{\mathrm{}}.$$
(14)


format: text

text:
# 1 S-Matrix Theory, Duality, and the Bootstrap

## 1 S-Matrix Theory, Duality, and the Bootstrap

In the late 1960s there were two parallel trends in particle physics. On the one hand, many hadron resonances were discovered, making it quite clear that hadrons are not elementary particles. In fact, they were found, to good approximation, to lie on linear parallel Regge trajectories, which supported the notion that they are composite. Moreover, high energy scattering data displayed Regge asymptotic behavior that could be explained by the extrapolation of the same Regge trajectories, as well as one with vacuum quantum numbers called the Pomeron. This set of developments was the focus of the S-Matrix Theory community of theorists. The intellectual leader of this community was Geoffrey Chew at UC Berkeley. One popular idea espoused by Chew and followers was “nuclear democracy” – that all hadrons can be regarded as being equally fundamental. A more specific idea was the “bootstrap”, that the forces arising from hadron exchanges are responsible for binding the hadrons, as composites of one another, in a more or less unique self-consistent manner.

The second major trend in the late 1960s grew out of the famous SLAC experiments on deep inelastic electron scattering. These gave clear evidence for point-like constituents (quarks and gluons) inside the proton. This led to Feynman’s “parton” model, which was also an active area of research in those days. I would submit that Susskind was philosophically more in tune with the parton world view than the S matrix one, though he clearly had a foot in each camp. My Berkeley training put me solidly in the S-matrix camp.

In the early 1970s, it became clear that QCD is the correct theory of strong interactions. It can be used for explicit computations in large momentum regimes where perturbation theory can be used, thanks to asymptotic freedom. Modern phenomenological studies in such regimes are directly descended from the older parton models. In other regimes, such as diffraction scattering, where perturbation theory is not applicable, S-matrix ideas are still used. So both major trends from the late 1960s have their descendents 30 years later.

String theory, which is the subject I want to focus on here, grew out of the S-Matrix approach to hadronic physics. The bootstrap idea got fleshed out in the late 1960s with the notion of a duality relating $`s`$-channel and $`t`$-channel processes that went by the name of “finite energy sum rules”  \- . Another influential development was the introduction of “duality diagrams”, which keep track of how quark quantum numbers flow in various processes . Susskind contributed papers on this topic . Later, duality diagrams would be reinterpreted as string world-sheets, with the quark lines defining boundaries. A related development that aroused considerable interest was the observation that the bootstrap idea requires a density of states that increases exponentially with mass, and that this implies the existence of a critical temperature, called the Hagedorn temperature  \- . It was generally supposed that this it is an ultimate temperature, though a phase transition was clearly another possibility.

## 2 The Dual Resonance Model

The duality program got a real shot in the arm in 1968 when Veneziano found a specific mathematical function that explicitly exhibits the features that people had been discussing in the abstract . The function, an Euler beta function, was proposed to give a good phenomenological description of the reaction $`\pi +\omega \pi +\pi `$ in the narrow resonance approximation. This was known to be a good approximation, because near linearity of Regge trajectories implies that the poles should be close to the real axis. A little later Lovelace and Shapiro proposed a similar formula to describe the reaction $`\pi +\pi \pi +\pi `$ . Chan and Paton explained how to incorporate “isospin” quantum numbers in accord with the Harari–Rosner rules . Also, within a matter of months Virasoro found an alternative formula with many of the same duality and Regge properties that required full $`s`$-$`t`$-$`u`$ symmetry . Later it would be understood that whereas Veneziano’s formula describes scattering of open-string ground states, Virasoro’s describes scattering of closed-string ground states.

In 1969 several groups independently discovered $`N`$-particle generalizations of the Veneziano four-particle amplitude  \- . The $`N`$-point generalization of Virasoro’s four-point amplitude was constructed by Shapiro . In short order Fubini and Veneziano, and also Bardakci and Mandelstam, showed that the Veneziano $`N`$-particle amplitudes could be consistently factorized in terms of a spectrum of single-particle states described by an infinite collection of harmonic oscillators  \- . This was a striking development, because it suggested that these formulas could be viewed as more than just approximate phenomenological descriptions of hadronic scattering. Rather, they could be regarded as the tree approximation to a full-fledged quantum theory. I don’t think that anyone had anticipated such a possibility one year earlier. It certainly came as a surprise to me.

One problem that was immediately apparent was that since the oscillators transform as Lorentz vectors, the time components would give rise to negative-norm ghost states. Everyone knew that such states would violate unitarity and causality. Virasoro came to the rescue by identifying an infinite set of subsidiary conditions, which plausibly could eliminate the negative-norm states from the spectrum . These subsidiary conditions are defined by a set of operators, which form the famous Virasoro algebra . The central term in the algebra was discovered by Joe Weis (unpublished). One price for eliminating ghosts in the way suggested by Virasoro was that the leading open-string Regge trajectory had to have unit intercept, and hence, in addition to a massless vector, it contributes a tachyonic ground state to the spectrum .

## 3 The String Idea

Once it was clear that we were dealing with a system with a rich spectrum of internal excitations, and not just a bunch of phenomenological formulas, it was natural to ask for a physical interpretation. The history of who did what and when is a little tricky to sort out. As best I can tell, the right answer – a one-dimensional extended object (or “string”) – was discovered independently by three people: Nambu, Susskind, and Nielsen. Nambu’s contribution to the 1969 conference held at Wayne State University apparently was first, but the conference was rather obscure , and the paper was not widely circulated. His subsequent paper, submitted to a symposium in Copenhagen, proposed that the string action is the area of the world-sheet in analogy with the proper length of the world-line of a point particle . These difficult to find papers have been reprinted in .

The string idea also appears quite clearly in Susskind papers, which were published in a refereed journal . Nielsen’s first paper on the subject was submitted to a 1970 conference in Kiev, though it was not published in the proceedings . However, shortly thereafter, he and Fairlie described their approach in a refereed journal . The Nielsen papers emphasize an analog electrostatic model in which one solves Laplace’s equation on a disk with sources on the boundary. This is just the string wave equation on a Euclideanized world sheet. The electrostatic analogy was also discussed by Shapiro for the Virasoro–Shapiro model taking the domain to be a sphere rather than a disk and with the sources attached to the interior of the surface . We recognize this to be the proper description of tree-level closed-string amplitudes.

A somewhat different approach treated the string world-sheet as some kind of sum or limit of complicated planar Feynman diagrams, sometimes referred to as fishnet diagrams. This approach, along with various related parton ideas, was pursued by Sakita and Virasoro  as well as by Susskind and Nielsen  \- . Note that these works are several years prior to ‘t Hooft’s famous paper on large-N gauge theory . They were not very specific, however, about which field theory should be used to form the fishnet diagrams or what limit was required to make contact with a string world sheet.

Though there was some progress in the intervening years , the string interpretation of the dual resonance model was not very influential in the development of the subject until the appearance of the 1973 paper by Goddard, Goldstone, Rebbi, and Thorn . It explained in detail how the string action could be quantized in light-cone gauge. Subsequently Mandelstam extended this approach to the interacting theory .

## 4 Loop Amplitudes

From the factorization of $`N`$-particle amplitudes one had learned the tree approximation spectrum and couplings. With this information in hand, it became possible to construct one-loop amplitudes. The first such attempt was made by Kikkawa, Sakita, and Virasoro . They did not have enough information in hand to do it completely right, but they pioneered many of the key ideas and pointed the way for their successors  \- .

Those who worked on this problem included a group at Princeton consisting of David Gross, André Neveu, Joël Scherk, and myself. Another group consisted of Lenny Susskind and various collaborators. Neveu and Scherk studied the divergence of the one-loop planar amplitude. They realized that by performing a Jacobi transformation of the theta functions in the integrand they could isolate the divergent piece and propose a fairly natural counterterm . Essentially the same thing was done independently by Susskind and Frye . The modern interpretation of these results is that viewed in a dual channel there is a closed string going into the vacuum. The divergence can be attributed to the tachyon in that channel. Its contribution is the piece they subtracted. This also explains why in a model without tachyons such divergences would not occur. The cancellation of the milder divergences due to dilaton tadpoles also became an important consideration in later years.

One of the important things discovered by the Princeton group , which was discovered independently by Frye and Susskind at about the same time , was that the nonplanar loop amplitude gives new and unexpected singularities. In both of these works the dimension of spacetime was assumed to be four, and the Virasoro subsidiary constraints were not implemented on the internal states circulating in the loop. As a result the singularities were found to be unitarity-violating branch points.

Lovelace observed that if one allows the spacetime dimension to be 26 and supposes that the subsidiary conditions imply that only transverse oscillators contribute, then instead of branch points the singularities would be poles . As we now know, these are are the closed-string poles in the nonplanar open-string loop. This calculation showed that unitarity requires that one choose the critical dimension and the intercept value for which the Virasoro conditions are satisfied. In those days these states were referred to as Pomeron-states rather than closed string states, since they necessarily carried vacuum quantum numbers. It was supposed that the unphysical Regge intercept value of two – implying the existence of a massless spin two particle as well as a tachyon – would somehow be lowered to the desired value of one in a more realistic model.

Shapiro carried out the first closed-string loop calculation . Specifically, he computed the one-loop (torus) amplitude.

## 5 Conclusion

The younger generation of theoretical physicists is aware of the wide range of important contributions that Lenny Susskind has made to our field in recent years. What they may be less aware of is his important role in the early history of string theory. I hope this presentation will help to publicize that work.


text:
# Directional coupling of optical signals by odd dark beams with mixed phase dislocations

## I Introduction

Dark spatial solitons (DSSs) formed in bulk self-defocusing Kerr-type nonlinear media (NLM) are able to guide co-propagating probe beams/pulses . The physical mechanism underlying is the intensity-dependent refractive-index change in a plane perpendicular to the propagation direction. Weak signal beams passing along these optically induced gradient waveguides  are subject to effective induced-phase modulation and are trapped. In photorefractive NLM  the situation is even more spectacular . Because of the wavelength-dependent material response, dark spatial solitons generated at low powers but at a photosensitive wavelength are able to guide and steer much more powerful ‘signal’ beams at non-photosensitive wavelengths. This is confirmed for photovoltaic , biased  and quasi-steady-state  photorefractive solitons and seems to hold also for incoherent dark solitons .

Bright signal beams/pulses can be deflected by both electronically  and optically induced prisms . When they are guided by DSSs, the same can be done by manipulating the transverse dynamics of the dark beams. The transverse velocity of an optical vortex soliton (OVS) has a radial and an angular component arising from the transverse phase and intensity gradients, respectively . Two practical ways to control the vortex rotation have their origin in the Guoy phase shift on both sides of the background beam waist  and in the interaction of ordered structures of OVSs  controlled by the topological charges. OVS steering is demonstrated by superposition of a weak incoherent background field . Operation of planar Y-junction splitters for signal beams is demonstrated in both Kerr-type  and photorefractive NLM  with pairs of grey DSSs born from even initial conditions. The possibility to branch a single input probe beam into ordered structures of sub-beams by quasi-two-dimensional DSSs is demonstrated numerically in Ref. . Other branching and steering schemes can be realized by employing the inherent dynamics of ring dark solitary waves , eventual NLM saturation  and/or anisotropy .

In this work we present numerical results on how one input optical data channel can be linked to a selected output data channel by means of an optically-induced steering waveguide. A high steering speed of the induced waveguide can be obtained by using odd dark beams (ODBs) of finite lengths  with a suitable choice of the mixed phase dislocation. Special attention is paid to ensure high energy efficiency of the directional coupler and short length of the interaction zone. Reconfiguration of the coupler is proposed to be done by changing the type of the phase dislocation reproduced by a multiple, active, single-voltage controlled computer-generated hologram .

## II Step-screw and edge-screw mixed phase dislocations

The mixed phase dislocations considered consist of a one-dimensional phase step of limited length, which ends, by necessity, with pairs of phase semi-spirals with opposite helicities. In this work two possible cases are considered. The first phase dislocation of this type, which will be denoted step-screw (SS), is described by the phase distribution

$$\mathrm{\Phi }_{\alpha ,\beta }^{SS}(x,y)=\mathrm{\Delta }\mathrm{\Phi }\left\{\frac{1\alpha }{2}\text{sgn}(y)\frac{\beta }{\pi }\mathrm{arctan}\left(\frac{\alpha y}{x+b\beta }\right)\right\},$$
(1)
where the parameters $`\alpha `$ and $`\beta `$ are defined as follows:

$$\alpha =\{\begin{array}{cc}0\hfill & \text{for }|x|b,\hfill \\ 1\text{ and }\beta =1\hfill & \text{for }x>b,\hfill \\ 1\text{ and }\beta =1\hfill & \text{for }xb.\hfill \end{array}$$
(2)
The phase distribution

$$\mathrm{\Phi }^{ES}(x,y)=\frac{\mathrm{\Delta }\mathrm{\Phi }}{2\pi }\left\{\mathrm{arctan}\left(\frac{y}{x+b}\right)\mathrm{arctan}\left(\frac{y}{xb}\right)\right\},$$
(3)
will be referred to as edge-screw (ES). In Eqs. (1)-(3) the quantity $`\mathrm{\Delta }\mathrm{\Phi }`$ stands for the magnitude of the step portion of the dislocation, $`2b`$ for its length, and $`x`$ and $`y`$ denote the transverse Cartesian coordinates parallel and perpendicular to the dislocation. Surface plots of the SS and ES phase dislocations are shown in Fig. 1a and 1b, respectively.

The formation of mixed phase dislocations was first identified by us in interferograms of decayed crossed one-dimensional DSSs in the presence of moderate saturation of the nonlinearity . From our present point of view the classification  of such formation generated via an instability as an ES or SS phase dislocation cannot be definitive. Both distributions shown in Fig. 1 are characterized by the same maximal phase difference $`\mathrm{\Delta }\mathrm{\Phi }`$ and dislocation length $`2b`$, but they differ in the rate $`\mathrm{\Phi }/\phi |_{|x|>b}`$ at which the phase changes azimuthally at their ends. In the step-screw case this change rate is twice as high as in the edge-screw case. In another work  we found experimentally and confirmed numerically that in the interval $`0.75\pi \mathrm{\Delta }\mathrm{\Phi }1.25\pi `$ the transverse steering velocity of an ODB of SS type is inversely proportional to $`\mathrm{\Delta }\mathrm{\Phi }`$. Because of the phase gradients across the edge portion of the dislocation the ES one starts ‘immediately’ steering. The step portion of the SS dislocation, however, is forced to steer by the outlying semi-helices. The process stabilizes after a propagation distance of the order of $`1L_{NL}`$ (see Figs. 2 and 3 in Ref. ), at a transverse velocity 2-3 times smaller than in the ES case. In view of this one more general statement holds: The transverse velocity of an ODB with mixed phase dislocation is inversely proportional to the absolute value of the azimuthal phase-change rate $`\mathrm{\Phi }/\phi |_{|x|>b}`$ outward its ends. Therefore the identification of the type of ODBs generated as a result of an instability is not very reliable.

## III Initial conditions for the model

The (2+1)-dimensional evolution of the ODBs in bulk homogeneous, isotropic and saturable NLM is described by the nonlinear Schrödinger equation

$$i\frac{E_D}{\zeta }+\frac{1}{2}\mathrm{\Delta }_{}E_D\frac{L_{Diff}}{L_{NL}}\frac{|E_D|^2E_D}{(1+s|E_D|^2)^\gamma }=0,$$
(4)
where $`\mathrm{\Delta }_{}=^2/\xi ^2+^2/\eta ^2`$. The transverse spatial coordinates are normalized to the initial dark beam width ($`\xi =x/a`$, $`\eta =y/a`$), and the nonlinear propagation path length $`\zeta `$ is expressed in units of Rayleigh diffraction lengths $`L_{Diff}=ka^2`$. Furthermore, $`L_{NL}=(kn_2I_0)^1`$ is the nonlinear length, $`k`$ is the wave number inside the NLM, and the background beam intensity at its entrance isnormalized to that needed to form an 1D DSS of width $`a`$ ($`I_0=I_{DSS}^{1D}`$). In order to maintain the correspondence to the conditions in our previous experimental work , where NLM with saturation was used, we adopted the refractive-index correction

$$\mathrm{\Delta }n=n_2|E_D|^2/(1+s|E_D|^2)^\gamma ,$$
(5)
with $`s=0.3`$ and $`\gamma =3`$. In a weak-signal approximation we took into account the bright-signal beam diffraction and the refractive-index changes it sees as induced by the ODB. In the NLM its evolution is described by

$$i\frac{E_B}{\zeta }+\frac{1}{2}\frac{\lambda _B}{\lambda _D}\mathrm{\Delta }_{}E_B\sigma \frac{L_{Diff}}{L_{NL}}\frac{|E_D|^2E_B}{(1+s|E_D|^2)^\gamma }=0.$$
(6)
The coupling coefficient $`\sigma `$ depends on the nature of the physical process of the optical nonlinearity (e.g. molecular orientation, electronic response of bound electrons, etc.) and on the experimental conditions (e.g. the polarizations of the waves). The numerical simulations presented in this work refer to the wavelength ratio $`\lambda _B/\lambda _D1`$ and $`\sigma =2`$. The slowly-varying electric-field amplitudes of the ODBs were assumed to be $`tanh`$-shaped and of the form

$$E_D(x,y,z=0)=\sqrt{I_0}B(r_{1,0}(x,y))\mathrm{tanh}\left[\frac{r_{\alpha ,\beta }(x,y)}{a}\right]e^{i\mathrm{\Phi }(x,y)},$$
(7)
where $`\mathrm{\Phi }(x,y)=\mathrm{\Phi }_{\alpha ,\beta }^{SS}(x,y)`$ for an SS, $`\mathrm{\Phi }(x,y)=\mathrm{\Phi }^{ES}(x,y)`$ for an ES phase dislocation, and $`r_{\alpha ,\beta }(x,y)=\sqrt{\alpha (x+\beta b)^2+y^2}`$ is the effective radial coordinate. In order to avoid any influence of the finite background beam of super-Gaussian form $`B(r_{1,0})=\mathrm{exp}\{(\sqrt{x^2+y^2/w^2})^{14}\}`$, its width was chosen to exceed the maximal ODB deflection at $`\zeta =10`$ more than 10 times. The signal beam is assumed to be $`sech`$-shaped and equal in width to the ODB width $`a`$. Its intensity is chosen hundred times weaker than $`I_{DSS}^{1D}`$, thus not disturbing the dark beams by cross-phase modulation. Equations (4) and (6) are solved numerically by modification of the beam propagation method over a $`1024\times 1024`$ grid.

Figure 2 shows the deflection of ODBs with edge-screw (ES) and step-screw (SS) phase dislocations vs. the nonlinear propagation path length $`\zeta =z/L_{NL}`$ for $`b/a=1.0`$ (solid), $`2.0`$ (dashed), and $`3.0`$ (dotted), normalized to the initial ODB width $`a`$. The difference between the transverse dynamics of the dark beams of different types is amazing. As already mentioned, it is due to the presence (in the ES case) or absence (in the SS case) of phase gradients across the one-dimensional portion of the dislocation at the entrance of the NLM. Besides by the type of the phase singularity, the transverse dynamics can be controlled in three possible ways : by chaging $`b/a`$ ratio, the phase step $`\mathrm{\Delta }\mathrm{\Phi }`$ and the background beam intensity. As it is seen in Fig. 2, changing the ODB length-to-width ratio gives an effective way to manipulate the ODB velocity. It is higher for the short and slower for the long dislocations. The long ones however could experience a transverse instability. More strongly pronounced is the increase of the velocity with decreasing the magnitude of the step portion of the phase dislocation $`\mathrm{\Delta }\mathrm{\Phi }`$ , but we refrained from exploiting this. The reason is that during the steering process both types of ODBs inevitably become grey (and asymptotically disappear), which tends to destroy the coupler modeled. The numerical calculations performed showed similar but slightly weaker coupling efficiencies when the transverse deflection is controlled by $`\mathrm{\Delta }\mathrm{\Phi }`$ than by $`b/a`$. All data in this work refer to $`\mathrm{\Delta }\mathrm{\Phi }(\zeta =0)=\pi `$. Characteristic of the evolution of the ES phase dislocation is that it disturbs the background beam more weakly than the SS dislocation. In order to emphasize this, the number of gray-scale levels in the frames shown in Fig 3 are intentionally reduced in the same way. The data refer to $`\zeta =z/L_{NL}=4`$ and $`b/a=1`$. The background beam intensity in this simulation is chosen to correspond to that of a 1D DSS ($`I_{DSS}^{1D}`$). The dashed line is intended to underline that the ES ODB is stronger deflected than the SS ODB.

The background beam intensity could be additionally used to controll the dynamics of the ODBs, even their deflection does not strongly depend on it . The background beam intensity, however, influences substantially the transverse beam profile and the guiding properties of the induced waveguide. Better spatial confinement could be expected at higher intensities. Numerical simulations were made to confirm this statement. Depending on the type of the dislocation (ES or SS) the highest encoupling efficiencies were achieved for levels $`(1.25÷1.75)I_{DSS}^{1D}`$. Higher intensities tend to destroy the induced waveguide due to the strong deformation of the dark beam profile under saturation of the nonlinearity.

The main task of the numerical simulations presented bellow is to find the proper parameter range for $`b/a`$ and beam intensity in order to reach maximal coupling efficiency in the analyzed schemes for directional coupling.

## IV One-directional coupler

As a first step in this work we modeled the coupling of an input signal beam (pulse) to a desired output channel by optically induced waveguides. Figure 4 shows the positions and the forms of the weak signal beams guided by the corresponding dark one after a nonlinear propagation path length $`\zeta =10`$. Channel $`0`$ is assumed to be supported by an OVS  with backround intensity $`1.5I_{DSS}^{1D}`$, channels $`1`$ and $`1`$ by ODBs with SS phase dislocations of $`b/a=1.5`$ at $`I=1.5I_{DSS}^{1D}`$, whereas channels $`2`$ and $`2`$ are considered to be induced by ODBs with ES phase dislocations with $`b/a=2.5`$ at $`I=1.25I_{DSS}^{1D}`$. The change of the ODB steering direction in the NLM (i.e. the switching between channels with positive and negative numbers) is controlled by reversing the gradients of the screw portions of the dislocations. First, real implementation of this directional switching is feasible if there are reliable means to change the phase distributions of the ODBs. The electrically-controllable, multiply-active, computer-generated holograms developed  are well suited for doing that. Because of the time needed to reconfigure the hologram, the relatively rare redirection of the data streams addressed to a desired channel will be preferable. Second, the energy efficiency of the guiding and switching should be as high as possible at the lowest possible cross-talk between the channels. The imaginary information channels considered to be coupled to the exit of the NLM are assumed to be rectangle-shaped with a width (height) 1.5 (3.0) times the input signal-beam FWHM. In Fig. 4 the respective positions along the steering direction (the arrow in Fig. 3) are sketched as rectangles.

The results for the directional coupling and cross-talk efficiencies obtained as ratios of the transmitted energy in a particular channel to the total one are summarized in Table I. The ODB with an SS phase profile is assumed to couple the input channel to channel No. $`1`$, whereas this with an ES dislocation redirects the input to channel $`2`$. The maximal efficiencies for channels $`0`$, $`\pm 1`$, and $`\pm 2`$ to be expected are $`74.7\%`$, $`70.7\%`$, and $`52.2\%`$, respectively. In most channels the cross-talk signal should remain well below $`14\%`$. The ODB steering is accompanied by the creation of a leading dark wave and a trailing bright peak on the background beam (see Fig. 3 in Ref. ) which are more strongly pronounced in the SS case. The bright trailing peak pushes the signal beam stronger. This qualitatively explains the relatively high efficiency of the coupling to channels $`+1`$ and $`1`$. Because of the additional self-defocusing of this peak it appears more effective to address channels $`+2`$ and $`2`$ by ES ODBs. If the out-lying channels are $`2.25`$ times larger than the input signal beam width $`a`$, in the ES case the efficiency in coupling the signal to them can be increased by $`3.6\%`$ to $`55.9\%`$. Due to the numerical discretization the overall accuracy in estimating the energies is within $`0.02\%`$ and the coupling efficiencies in the respective channels are correct to within $`\pm 0.04\%`$.

We should note here that because of the saturation of the nonlinearity the width of the induced waveguide is larger than in the case of pure Kerr nonlinearity. If a Kerr medium is considered (i.e. $`s=0`$) the coupling efficiencies increase by $`10\%`$, whereas the cross-talk coefficients become slightly weaker.

## V Two-directional coupler

For addressing a desired channel from an imaginary array of information channels ($`3\times 3`$ in this work) we used the squared-shaped channel configuration shown in Fig. 5. The data presented below refer to $`\zeta =10`$ and channel width and height $`1.5`$ times the input signal beam FWHM. Best spatial confinement of the signal beams is achieved at $`1.5`$ times higher background beam intensity than that needed to form a 1D DSS of width $`a`$. Besides being controlled by the type of the phase dislocation (screw, ES, SS), the ODB transverse velocity and deflection at the exit of the NLM are adjusted by the initial ODB length-to-width ratio $`2b/a`$ ($`b/a=1.5`$ in the SS case; and $`b/a=3.4`$ in the ES case; both at $`\mathrm{\Delta }\mathrm{\Phi }=\pi `$). Let us denote the channels as $`C_{IJ}`$, where $`I`$ and $`J`$ are the respective row and column indices. The central output channel $`C_{22}`$, which is on axis with respect to the input one, is addressed by an OVS , the side channels $`C_{21}`$, $`C_{23}`$, $`C_{12}`$, and $`C_{32}`$ by ODBs with appropriately $`\pm 90^{}`$ rotated SS phase dislocations, whereas signal beam coupling to each of the diagonal channels $`C_{11}`$, $`C_{13}`$, $`C_{31}`$, and $`C_{33}`$ is performed by ODBs with ES phase dislocations rotated $`\pm 45^{}`$. The calculated coupling efficiencies are summarized in Table II. The maximal efficiencies for the direct channel $`C_{22}`$, the neighboring channel $`C_{23}`$, and the diagonal channel $`C_{11}`$ to be expected are $`61.6\%`$, $`49.6\%`$, and $`47.1\%`$ respectively, at cross-talk signals remaining well below $`12\%`$. In Fig. 5 we composed grey-scale frames of the signal beams as located on the imaginary array of information channels considered. Coupling the input channel to information channel $`C_{11}`$ (by an ES ODB) should result in undesired but equal cross-signals in channels $`C_{IJ}`$ and $`C_{JI}`$ ($`IJ`$). If the input signal is coupled to channel $`C_{23}`$ the same should hold for channels $`C_{1J}`$ and $`C_{3J}`$. The data obtained indicate that the discretization in our simulations leads to a maximal inaccuracy in estimating the coupling efficiencies that is within $`\pm 0.04\%`$. If pure Kerr nonlinearity is considered the coupling efficiencies are $`1\%`$ to $`8\%`$ higher.

## VI Conclusion

The numerical data on the nonlinear evolution of mixed step-screw and edge-screw phase dislocations under equivalent conditions shows an amazing difference between the transverse steering dynamics of SS and ES ODBs. Such odd dark beams are able to induce steering conduits for weak information beams (pulses) in bulk nonlinear media. Due to this two different schemes for directional coupling are proposed. The relatively high coupling and moderate cross-talk efficiencies in both of them provide a reasonable basis for further optimization and may open the way to constructing parallel reconfigurable all-optical logical elements.

###### Acknowledgements.

D.N. thanks Vrije Universiteit Amsterdam for the opportunity to complete part of the computational work at the facilities of the Laser Center, VU, Amsterdam. A.D. is grateful to the Alexander von Humboldt Foundation (Germany) for the award of a fellowship and the opportunity to work in the stimulating atmosphere of the Max-Planck-Institut für Quantenoptik (Garching, Germany).


text:
# To Probe a Black Hole

## 1 Introduction

Sagittarius A\* is the largest black hole as seen from Earth, with $`2R_S/D=15\mu \mathrm{as}`$, where $`R_S`$ is the Schwarzschild radius and $`D`$ is the distance to the Galactic Center (Genzel et al 2000, Ghez et al 1998). Several other nearby supermassive black holes have angular diameters $`10\mu \mathrm{as}`$ (Richstone et al 1998).

$`15\mu \mathrm{as}`$ is the angular resolution of a 375 GHz VLBI, and the interstellar scattering size towards Sgr A\* at 375 GHz is only $`6\mu \mathrm{as}`$. Sgr A\* might have already been detected by a 215 GHz VLBI (Krichbaum et al 1998). Radio-astronomers believe that future observations will resolve the black hole (Doeleman & Krichbaum 1999, Falcke et al 1999). The high-frequency VLBI radio map can show the light capture cross section (the angular diameter is $`38\mu \mathrm{as}`$ 111For simplicity, here and in the rest of the paper, we assume that the black hole is not rotating. For concreteness, we discuss Sgr A\*.), the last stable orbit ($`45\mu \mathrm{as}`$), the innermost jet formation region, or the combination of the above.

When the $`40\mu \mathrm{as}`$ structure is detected, one will want to use the radio map to probe the black hole metric. The standard view is that this is impossible. What we see is determined by both the metric and the emissivity. The emissivity depends on the plasma dynamics. The dynamics of a magnetized collisionless plasma in a black hole metric will not be understood in a foreseeable future. Even in Newtonian approximation, our understanding of the low-luminosity accreting black holes is rudimentary (Shvartsman 1971, Gruzinov 1997, 1999, Quataert & Gruzinov 1999, 2000 (a,b)).

Nevertheless, an idealized high-frequency VLBI can probe the black hole metric. This is possible if the emissivity is variable in space and time. Lensing of this variable emission should be detectable, and can be used to probe the metric.

In §2 we discuss one of many possible ways to detect lensing, and therefore to probe the black hole. Only real data can tell which method would work best. The aim of this paper is to show that high-frequency VLBI can probe black holes in a plasma-dynamics-independent way – in principle.

## 2 Lensed flares

The size of Sgr A is known to be small. The observed size is scattering-dominated at lower frequencies, up to about 100 GHz. From 215 GHz VLBI, the radius of the emitting region is thought to be in the range 3 – 13 $`R_S`$. The flux probably peaks at about 1000 GHz (Zylka et al 1992, Serabyn et al 1997). We will assume that most of the high-frequency radio emission is from a few $`R_S`$ region around the hole. In what follows we make some further assumptions about the emissivity. These assumptions are only qualitative, the suggested diagnostic does not use any quantitative model of the plasma flow.

Define a radio flare as a sudden brightening of a small region. Here sudden means $`R_S/c0.5\mathrm{min}`$ and small means $`R_S/D10\mu \mathrm{as}`$. We will assume that VLBI has a comparable angular and time resolution 222Good uv-coverage is a plus but not a must..

Every flare is lensed by the hole an infinite number of times, but the flux even in the second image is likely to be very small if the source is far from the hole. The situation is different for sources close to the hole. At $`4R_S`$ say, the probability that the magnification of the second image is greater than 0.2 is about 10% (Appendix A).

If flares do occur, lensed double flares will be detected. If the background emission is detected in the same experiment, one can fix a position (or an assumed position) of the hole in the sky. Given relative positions of the two images and the hole, one can calculate the theoretical time lag (Appendix A). We probe the metric by comparing the theoretical and the observed time lags. For example, a source at $`R=5R_S`$ and $`\alpha =0.2`$ (fig. 1) will be seen as a double (magnifications $`2`$ and $`1`$). The brighter image is 32$`\mu `$as from the hole, the second image is 28$`\mu `$as from the hole. The second image appears 50s later.

###### Acknowledgements.

This work was supported by the W. M. Keck Foundation and NSF PHY-9513835. I thank Eliot Quataert, John Bahcall, and Avi Loeb for useful discussions.

## Appendix A Lensing in the Schwarzschild metric

We use dimensionless units, $`R_S=c=1`$. In Schwarzschild spherical coordinates $`(r,\theta ,\varphi )`$, the observer is at $`(D,0,0)`$, where $`D1`$. The source is at $`(R,\pi \alpha ,0)`$, as shown in figure 1. Only the two brightest images are considered.

Light propagation is given by (eg. Shapiro & Teukolsky 1983)

$$\dot{r}=b(1r^1)\sqrt{b^2r^2+r^3},$$
(A1)
$$\dot{\psi }=b(1r^1)r^2.$$
(A2)
Where the overdot is the time-at-infinity derivative, $`b`$ is the impact parameter, and $`\psi `$ is the polar angle in the plane of motion. The impact parameters of the two brightest images are given by (fig.1.):

$$\pi \pm \alpha =\left(_0^{x_2}+_{R^1}^{x_2}\right)\frac{dx}{\sqrt{b^2x^2+x^3}},$$
(A3)
where $`x_2`$ is the intermediate root of the polynomial $`b^2x^2+x^3`$. The minus sign in (A3) corresponds to the brighter image. The propagation time to infinity is

$$t=\left(_0^{x_2}+_{R^1}^{x_2}\right)\frac{dx}{b(1x)x^2\sqrt{b^2x^2+x^3}}.$$
(A4)
The integral diverges, but the time lag $`t=t_2t_1`$ is finite. The magnification is given by the ratio of the solid angles of emitted and lensed rays,

$$\mu =\frac{d(\mathrm{cos}\beta )}{d(\mathrm{cos}\alpha )},$$
(A5)
where the emission angle $`\beta `$, as measured by an observer at rest at $`R`$, is given by (Shapiro & Teukolsky 1983)

$$\mathrm{sin}\beta =bR^1\sqrt{1R^1}.$$
(A6)


text:
# I Introduction

## I Introduction

The occurence of oscillations and chaos in neural networks  is receiving a growing interest. From a neurophysiological point of view, chaotic behavior has been observed in brain dynamics  and dynamics reduction mechanisms have been proposed to explain experimental findings concerning cognitive processes in the brain . The study of complex dynamics in simulated neural networks started with the paper by Sompolinsky and coworkers , where it was pointed out that asymmetric synapses may lead to chaos; subsequently chaotic neural networks have been investigated in a number of papers, all focusing on the role played by random asymmetric synapses .

In the last years several examples of chaotic dynamics arising from the choice of a nonmonotonic transfer function (i.e. the function that gives the state of the neuron as a function of the postsynaptic potential) have been presented . The physiological justification of a nonmonotonic transfer function is ascribed to the fatigue of neurons after being exposed to a large post-synaptic input . In  a strongly diluted neural network with nonmonotonic binary neurons and adapting synapses was studied and dynamics reduction was observed for low connectivity. The interplay between microscopic and macroscopic dynamics in this class of networks was analyzed in : the phenomenon of damage spreading showed that at the microscopic level the network’s motion is always to be considered chaotic, whereas from a macroscopic, mean-field point of view a rich variety of behaviors can occur: fixed-point, periodic attractors, and strange attractors.

Considering neural networks in the more general frame of input-output devices, neural nets with nonmonotonic transfer functions have been employed in applications such as generation of robust chaos  and clustering of data .

The purpose of this paper is to extend the above mentioned results to the case of mean-field networks of continuous neurons and synapses. Moreover the influence of noisy external stimuli as a mechanism for controlling the chaotic dynamics  in this kind of neural networks is analyzed.

The model consists of $`N`$ continuous neurons $`S_i[1,+1],i=1,\mathrm{},N`$; for each neuron $`S_i`$, $`K`$ input sites $`j_1(i),\mathrm{},j_K(i)`$ are randomly chosen among the $`N`$ sites, and $`NK`$ real valued synaptic couplings $`J_{ij}(t)`$ are introduced. The local field acting on neuron $`S_i`$ at time $`t`$ is $`h_i(t)=J_{ij}(t)S_j(t)`$, the sum taken over the $`K`$ random input neurons. Neurons and synapses evolve according to a parallel discrete time deterministic dynamics, the updating rule for neurons being

$$S_i(t+1)=f(h_i(t)),$$
(1)
where $`f`$ is a nonmonotonic transfer function whose form is the following:

$$f(x)=\{\begin{array}{cc}x/\theta & \text{if }|x|<\theta \hfill \\ \text{sign}(x)& \text{if }\theta <|x|<c\theta \hfill \\ 0& \text{if }|x|>c\theta .\hfill \end{array}$$
(2)
For low values of local field, $`f`$ is linear with gain $`\frac{1}{\theta }`$, for intermediate values it reduces to the $`sign`$ function; eventually, for large postsynaptic inputs, it pulls the neuron into the quiescent state ($`S=0`$), thus modeling the fatigue phenomenon. The nonmonotonicity of $`f`$ is controlled by the threshold $`\theta `$, the constant $`c`$ (we fix $`c=2`$ in the following) determining the width of the range of $`h`$ values which saturate neurons.

A Hebbian evolution law for synapses is assumed as follows:

$$J_{ij}(t+1)=(1A)J_{ij}(t)+AS_i(t)S_j(t),$$
(3)
where $`A`$ is the ’learning rate’. In order to preserve network’s functioning, the value of $`A`$ must be sufficiently small so that the learning speed is not too large (neural networks with adapting synapses are discussed in  and several papers cited therein). Obviously, the fixed-synapses situation can be recovered by fixing $`A`$ to zero.

In the limit of large $`N`$, while keeping fixed the connectivity $`K`$, one obtains a strongly diluted network  , whose dynamics is described by mean-field equations (see ). In other words, one can treat neurons and synapses as independent stochastic variables and assume a Gaussian distribution for the local field acting on neurons. The following macroscopic parameters are introduced to describe the neuronic configuration: the overlap with an arbitrary pattern $`\{\xi \}`$ (we choose $`\{\xi _i=1\}`$ for simplicity) $`m(t)=S_i(t)`$ and the neuronic activity $`q(t)=S_i^2(t).`$ We remark that the suppression of the site index is possible because all averages are site independent. The synaptic configuration is described by the average strength $`J(t)=J_{ij}(t)`$ and the activity $`W(t)=J_{ij}^2(t)`$. Mean and variance of the local field are then given respectively by:

$$\begin{array}{cc}\mu (t)=Km(t)J(t),& \\ & \\ \sigma ^2(t)=K(W(t)q(t)J^2(t)m^2(t)).& \end{array}$$
(4)
Simple calculations provide the flow equations for the four macroscopic parameters of the network:

$$m(t+1)=𝑑h\frac{1}{\sqrt{2\pi \sigma ^2}}\text{exp}\left[(h\mu )^2/2\sigma ^2\right]f(h)$$
(5)
$$q(t+1)=𝑑h\frac{1}{\sqrt{2\pi \sigma ^2}}\text{exp}\left[(h\mu )^2/2\sigma ^2\right]f^2(h)$$
(6)
$$J(t+1)=(1A)J(t)+Am^2(t)$$
(7)
$$W(t+1)=(1A)^2W(t)+2A(1A)J(t)m^2(t)+A^2q^2(t)$$
(8)
In what follows we shall also study the mean field behavior of the system when an external noisy stimulus is presented to the network. We shall assume that stimuli are independent Gaussian variables with zero mean and time-dependent variance $`I(t)`$. The only effect of these stimuli is then to increase the fluctuations of local field around its mean, indeed the variance of $`h(t)`$, in presence of external stimuli, reads:

$$\sigma ^2(t)=K\left[W(t)q(t)J^2(t)m^2(t)\right]+I(t).$$
(9)
In the next section some properties of flow equations, in the absence of external stimuli, are shown, while in Section 3 the influence of external periodic stimuli acting on the network is considered. Section 4 summarizes our conclusions.

## II Dynamics in the absence of external stimuli

Let us start by considering the system’s motion in the absence of external stimuli and perform the analysis which has been applied in  in the case of binary neurons.

Firstly we consider $`A=0`$, i.e. fixed synapses: eqs. (7) and (8) can be dropped in this case and the time evolution of the system is described by eqs. (5) and (6) for the parameters $`m`$ and $`q`$ which characterize the neuronic configuration. The structure of these two equations is similar to the two-dimensional map studied in  and we find the same qualitative behavior. The macroscopic dynamics is represented in the bifurcation diagram of Fig. 1: depending on the values of parameters the macroscopic trajectory can be chaotic, periodic or converge to a fixed point. We verified that the system becomes chaotic via the Feigenbaum scenario, the bifurcation mechanism being period doubling. In Fig. 2 (left) we depict the strange attractor corresponding to $`K=15`$, $`J=0.8`$, $`W=0.9`$ and $`\theta =3`$: an unstable fixed point lies on the attractor. We discuss now the hyperbolicity of the attractor; we remark that in the hyperbolic case many interesting properties about the structure and dynamics of chaos have been demonstrated . In Fig. 2 (right), the stable manifold of the fixed point, i.e. the set of points $`\{m,q\}`$ such that the orbit starting from $`\{m,q\}`$ approaches the fixed point, is depicted. Since the stable manifold shows near tangencies with the attractor, we conclude that the attractor is not hyperbolic: a similar conclusion was given in  with respect to binary neurons, therefore nonhyperbolicity is a general feature of the macroscopic chaos shown by this class of neural networks.

Lyapunov exponents give a means of characterizing the stretching and folding characteristics of attractors. We evaluate them by the method described in ; for the case depicted in Fig. 2 we obtain $`\lambda _1=0.65`$ and $`\lambda _2=3.23`$. Having checked that always at most one Lyapunov exponent is positive, we conclude that the macroscopic chaos provided by this system is fragile , i.e. a slight change of the parameters can typically destroy chaos so that a stable periodic orbit sets in. Concerning the fractal dimension of the attractor, we follow the method proposed by Grassberger and Procaccia . A finite set of time-delayed vectors is constructed from a sample of $`N`$ points on the attractor, and the correlation integral $`C`$, at scale $`r`$, is estimated by counting the number of pairs of vectors whose distance does not exceed $`r`$. The scaling law $`C(r)r^\nu `$ allows to measure the exponent $`\nu `$ which is related to the Hausdorff dimension $`D`$ by $`\nu D`$, and can be considered itself a measure of the strangeness of the attractor. In Fig.3 the log-log plot of the correlation integral versus the scale $`r`$ is reported for the case in Fig. 2: we evaluate $`\nu =1.07`$.

We now turn to consider adapting synapses, fixing $`A`$ small and greater than zero. We confirm the results obtained in : for low connectivity chaos is removed in the stationary regime. In Fig. 4 the bifurcation diagram is depicted for $`K=8`$: only periodic attractors arise. On the other hand, we find that chaos is not completely removed for high connectivity, in agreement with .

The discussion presented in this section may be summarized concluding that the features of macroscopic chaos displayed by diluted neural networks with nonmonotonic transfer function do not depend on the nature of the neurons, since we find the same qualitative behavior whether the neurons are binary or real valued. In particular in  the role of the fluctuations of local field in determining the dynamics reduction was outlined: since the presence of noisy external stimuli leads to increasing fluctuations of $`h(t)`$, as already discussed in the previous section, we expect that a proper stimulus might control the chaotic behavior of these systems. The next section deals with this issue.

## III External stimuli

The classical strategy for chaos control is the well known Ott-Grebogi-Yorke method , where unstable periodic orbits are stabilized by means of small external perturbations to system parameters. Another interesting approach has been proposed in  and is based on the application of periodic feedbacks acting on system variables (instead of system parameters). The application of this method to small size neural systems was described in , where the stabilization mechanism was presented as a technical tool for chaos control, without a physical motivation to support the choice of the form of the perturbation. In the following we provide a physical justification for chaos-control strategies like the one proposed in . We show that in the class of neural networks here considered a mechanism for chaos control exists which has a natural physical interpretation, i.e. periodic perturbations by noisy external stimuli acting on the network.

In the following we limit our discussion to the fixed synapses case. We fix the parameters of the network so that, in the absence of stimuli (autonomous system), the trajectory of the system would be chaotic. A noisy external stimulus of variance $`I`$ is applied to the network every $`p`$ time steps, i.e. its time evolution is assumed to be:

$$I(t)=I\delta _{t,p},$$
(10)
where $`\delta _{t,p}=1`$ when $`t`$ is a multiple of $`p`$ and zero otherwise. The stationary regime of the perturbed system is described in the bifurcation diagrams of Fig. 5 which clearly illustrate how the chaotic dynamics is controlled by the noisy stimulus: the cases $`p=2`$ (a), $`p=3`$ (b) and $`p=4`$ (c) are shown. For low values of $`I`$ the network remains chaotic; as $`I`$ is increased, the chaotic region intercalates with windows of periodicity. By successive bifurcations the system enters in a wide window characterized by a periodic behavior with the same period as the external stimulus. As $`I`$ is further increased, a continuous transition from the $`p`$-periodic behavior to the trivial fixed point $`m=0`$ is observed. We have thus shown that there exists a wide range of $`I`$-values which force the system to attain a periodic orbit with the same period as the stimulus.

Another interesting situation corresponds to periodic noisy stimuli with large duration, like the one whose time evolution is depicted in the lower part of Fig. 6. Also in this case the stimulus forces the system to attain a periodic orbit with the same period of the stimulus (upper part of Fig. 6). Interestingly we find that during the time intervals corresponding to $`I>0`$ the network’s variables $`\{m,q\}`$ are almost constant with values very close to those corresponding to the unstable fixed point lying on the attractor of the autonomous system: external periodic perturbations seem, in this case, to stabilize the unstable fixed point of the network. We obtain similar results varying the parameters of the autonomuous systems and the width of pulses in the waveform of $`I(t)`$.

## IV Conclusions

In this paper we have studied diluted neural networks with continuous neurons and nonmonotonic transfer function. We have shown that in passing from the case of binary neurons (treated in previous papers) to continuous neurons, the dynamics is described by the same qualitative scenario. In particular, the adapting system removes chaos, in the stationary regime, for low connectivity. We also outlined a mechanism for controlling chaos in these neural systems, consisting in submitting the network to a noisy stimulus whose variance is periodic. For a proper amount of external perturbation the system is forced to behave periodically with the same period as the stimulus.

## Ackowledgements

The authors thank L. Angelini, C. Marangi and M. Pellicoro for useful discussions.

Figure Captions

Figure 1: Bifurcation diagram for the system with fixed synapses, while keeping fixed $`K=15`$, $`J=0.8`$, $`W=0.9`$.

Figure 2: (Left) The attractor of the system with fixed synapses corresponding to $`K=15`$, $`J=0.8`$, $`W=0.9`$ and $`\theta =3`$. (Right) The stable manifold of the unstable fixed point lying on the attractor (put in evidence by a circle centered on it) is depicted togheter with the attractor.

Figure 3: Log-log plot of the correlation integral (see  for details) versus the scale $`r`$ for the attractor depicted in Fig. 2. The slope of the curve provides the Hausdorff dimension.

Figure 4: Bifurcation diagram for the adapting system. The parameters are: $`K=8`$, $`J=0.8`$, $`W=0.9`$, $`A=10^4`$.

Figure 5: Bifurcation diagrams for the system perturbed by a noisy stimulus with period $`p`$ and variance $`I`$. The parameters are $`K=15`$, $`J=0.8`$, $`W=0.9`$ and (a) $`p=2`$, (b) $`p=3`$, (c) $`p=4`$.

Figure 6: The time evolution of the system with fixed synapses in presence of noisy periodic stimuli. The parameters are $`K=15`$, $`J=0.8`$, $`W=0.9`$. (Top) Time evolution of the overlap $`m(t)`$. (Middle) Time evolution of the neuronic activity $`q(t)`$. (Bottom) The waveform representing the time evolution of the variance of the external stimulus: each pulse is 15-time steps long.