Boris Kartascheff
email: kartascheff@fkp.tu-darmstadt.de

Christian Guill
email: guill@fkp.tu-darmstadt.de

Barbara Drossel
email: drossel@fkp.tu-darmstadt.de

Positive complexity-stability relations in food web models without foraging adaptation

Keywords: food web structure, niche model, vulnerability distribution, foraging efforts, functional response

May's local stability analysis of random food web models showed that increasing network complexity leads to decreasing stability, a result that is contradictory to earlier empirical findings. Since this seminal work, research of complexity-stability relations became one of the most challenging issues in theoretical ecology. We investigate conditions for positive complexity-stability relations in the niche, cascade, nested hierarchy, and random models by evaluating the network robustness, i.e. the fraction of surviving species after population dynamics. We find that positive relations between robustness and complexity can be obtained when resources are large, Holling II functional response is used and interaction strengths are weighted with the number of prey species, in order to take foraging efforts into account. In order to obtain these results, no foraging dynamics needs to be included. However, the niche model does not show positive complexity-stability relations under these conditions. By comparing to empirical food-web data, we show that the niche model has unrealistic distributions of predator numbers. When this distribution is randomized, positive complexity-stability relations can be found also in the niche model.

A c c e p t e d m a n u s c r i p t

Positive complexity-stability relations in food web models without foraging adaptation

Boris Kartascheff 1 , Christian Guill 2 , and Barbara Drossel 3 Institute of Solid-state Physics, Darmstadt University of Technology, Hochschulstraße 6 D-64289 Darmstadt, Germany

A c c e p t e d m a n u s c r i p t 1 Introduction

Diverse ecosystems such as tropical rainforests are apparently far more stable than monocultures in agriculture, which are vulnerable to pest invasions. Based on this and other observations, the traditional view of ecologists was formed that more complex ecosystems are more stable, see e.g. the publications by [START_REF] Odum | Fundamentals of ecology[END_REF], [START_REF] Macarthur | Fluctuations of animal populations and a measure of community stability[END_REF], [START_REF] Elton | Ecology of Invasions by Animals and Plants[END_REF], and others in the second half of the 20th century. This view was overturned when [START_REF] May | Will a large complex system be stable?[END_REF] showed that complex community models are inherently unstable. Based on random matrix theory, he analytically found conditions for a dynamical equilibrium state to be stable: Let α be the average interaction strength. If the system has Lotka-Volterra dynamics, the interaction strength between two species is proportional to the coupling coefficient in the dynamical equations. Let S be the number of species in the system and C the connectance, which is defined as C = L/S 2 , with L being the number of nonzero connections. According to May, a network is "almost certainly stable" if

α < (S • C) -1/2 (1)
and "almost certainly unstable" if

α > (S • C) -1/2 .
(2)

May's findings had a remarkable impact on theoretical ecology, giving birth to the "complexity-stability debate" [START_REF] Mccann | The diversity-stability debate[END_REF][START_REF] Garcia-Domingo | Food-web complexity emerging from ecological dynamics on adaptive networks[END_REF]. This debate tries to reconcile the fact that empirical foodwebs span a wide range of values of S and C (see, e.g., the examples listed in [START_REF] Dunne | Food web structure and network theory: the role of connectance and size[END_REF]) and thus appear not to be limited by an upper bound on the product S • C, with the results obtained by May.

May's analysis never stood uncriticized: First, the structure of real food webs is far from random, as is well known from food-web topology studies [START_REF] Williams | Simple rules yield complex food webs[END_REF][START_REF] Dunne | Food web structure and network theory: the role of connectance and size[END_REF]. Indeed, when link strengths and topology are carefully selected either by copying them from real food webs [START_REF] Yodzis | The stability of real ecosystems[END_REF] or by assembling a web by repeated addition of new species [START_REF] Law | Self-Assembling Food Webs: A Global Viewpoint of Coexistence of Species in Lotka-Volterra Communities[END_REF][START_REF] Morton | Regional species pools and the assembly of local ecological communities[END_REF], one can generate stable complex webs. Second, food webs do not need to be at an equilibrium and need not return to the previous state after a small perturbation in order to show persistence of all species. For this reason, several authors use a stability measure which evaluates the robustness (R) of entire communities (i.e. the fraction of species that survive after population dynamics). Such robustness investigations are usually performed with more realistic food web structures than that of random webs, and they incorporate nonlinear population dynamics. In order to search for conditions under which robustness can remain large with increasing S • C, these investigations evaluate robustness either as function of S at fixed C, or as function of C at fixed S. However, even when the mentioned more realistic systems are used, the data typically show a negative complexity-stability relation, i.e., robustness decreases when either S or C is increased (Kondoh, 2003a;[START_REF] Garcia-Domingo | Food-web complexity emerging from ecological dynamics on adaptive networks[END_REF][START_REF] Kondoh | Does foraging adaptation create the positive complexity-stability relationship in realistic food-web structure?[END_REF][START_REF] Uchida | Relation between complexity and stability in food webs with adaptive behavior[END_REF].

Recently, it was found that including foraging adaptation can lead to positive complexitystability relations when random or cascade topologies are used (Kondoh, 2003a). Foraging adaptation is the active adaptation of prey preferences in order to maximize total

A c c e p t e d m a n u s c r i p t

food uptake, in contrast to passive switching between alternative prey species that simply arises due to density dependent interaction rates [START_REF] Gentleman | Functional responses for zooplankton feeding on multiple resources: a review of assumptions and biological dynamics[END_REF]. Adaptive foraging is implemented by allowing predators to assign more search time to prey that is more abundant. Additional dynamical equations are included in the models for the foraging efforts, which increase (decrease) when the gain obtained per unit effort from a given prey is larger (smaller) than the average gain obtained from all prey, thus maximizing the food intake of a species. A positive complexity-stability relation was found for increasing C (but not for increasing S) in the random and cascade models, but not in the niche model [START_REF] Brose | Comment on "Foraging adaptation and the relationship between food-web complexity and stability[END_REF]. The reason is that the number of basal species decreases in the niche model when C is increased for fixed S. This is in contrast to other models, where the number of basal species increases with increasing C. When the niche model is modified such that the number of basal species is kept constant [START_REF] Kondoh | Does foraging adaptation create the positive complexity-stability relationship in realistic food-web structure?[END_REF][START_REF] Uchida | Relation between complexity and stability in food webs with adaptive behavior[END_REF], a positive complexity-stability relation is also found in the niche model with foraging adaptation. When S is increased for fixed C, robustness decreases in all these models. It is the purpose of this paper to show that even without foraging adaptation positive complexity-stability relations can be obtained. We simulate foodwebs that are created by the rules of the random, cascade [START_REF] Cohen | A stochastic theory of community food webs. Models and aggregated data[END_REF], niche [START_REF] Williams | Simple rules yield complex food webs[END_REF] and nested hierarchy model [START_REF] Cattin | Phylogenetic constraints and adaptation explain food-web structure[END_REF]. The advantage of these relatively simple models is that they need as an input only the two parameters S and C, which are directly related to the complexity of the networks. Furthermore, since these models provide a direct rule for constructing a network, they can be used to produce networks with distinctive topological features. This allows us to investigate the effect of various topological features on robustness under population dynamics. More advanced models that aim at reproducing empirical food web topologies are available, such as the speciation model [START_REF] Rossberg | An explanatory model for food-web structure and evolution[END_REF] or the matching model [START_REF] Rossberg | Food webs: Experts consuming families of experts[END_REF]. While these models reproduce empirical data even better than the already very good niche or nested hierarchy model, they are not suitable for our analyzes mainly for two reasons: First, they need as an input up to eight parameters that are used to fit the model to empirical data and that cannot be directly related to topological features. Second, they construct food webs by an evolutionary algorithm. This has the merit of mimicking the dynamical processes that lead to empirically observable food web structures, but makes it difficult to deduce the interplay between topological characteristics of the model networks and their robustness. We obtain a positive relation for increasing C under the conditions that the topology is of random, cascade or nested hierarchy type, that resource biomass is sufficiently large, that the functional response is of Holling type II [START_REF] Holling | Some characteristics of simple types of predation and parasitism[END_REF], and that foraging efforts are taken into account by weighting the interaction strengths with the number of prey species. These foraging efforts are implemented by dividing the encounter rate of a predator with one of its prey by the number of prey species of this predator, such that generalist species, which have many prey species, can only spend part of their available time to search for one particular prey. For the niche model, there is no positive complexity-stability relation under these conditions, and we found that this is due to unrealistic distributions of the number of predator species per prey species (vulnerability distributions) in the niche model. We show that the vulnerability distributions of real food webs are different from those of the niche model. When the vulnerability distributions are randomized in the model, we find positive complexity-stability relations.

A c c e p t e d m a n u s c r i p t

For the nested-hierarchy model we find positive complexity-stability relations when C is increased, but not when S is increased. The latter trend can be inverted by randomizing vulnerability distributions, but since in the nested hierarchy model these distributions fit empirical data quite well, the randomization makes the model less realistic and has therefore no biological meaning.

Model

Food-web structure

Model food webs are constructed by assigning links at random within the constraints given by the rules of the cascade [START_REF] Cohen | A stochastic theory of community food webs. Models and aggregated data[END_REF], niche [START_REF] Williams | Simple rules yield complex food webs[END_REF], nested hierarchy [START_REF] Cattin | Phylogenetic constraints and adaptation explain food-web structure[END_REF] , and random models. The parameters for each of these models are the number of species S and the connectance C. They determine the complexity of a network. Two resources with constant and equal size Res are assigned to each model web. They serve as an inexhaustible pool of nutrients, and since they do not depend on the population dynamics of the trophic species, they can be regarded as an external source of energy of the food web. In the niche and nested hierarchy models, the resources were assigned the smallest niche values. In niche model food webs, species with a feeding range that does not cover the resources are nevertheless assigned a link to the resources if their feeding range contains no species at all, in accordance with the original algorithm by Williams and Martinez [START_REF] Williams | Simple rules yield complex food webs[END_REF]. In the random and cascade models, links to the resources are assigned in the same way as all other feeding links.

Including the external resources as special nodes in the model food webs is a simple way of obtaining a variable number of basal species (i.e., trophic species that feed from at least one of the resources) according to the rules of the corresponding model algorithm rather than fixing it to an arbitrary number. When investigating the niche model, some authors use a modification where the number of basal species (species that feed from the resource) is kept constant [START_REF] Kondoh | Does foraging adaptation create the positive complexity-stability relationship in realistic food-web structure?[END_REF][START_REF] Brose | Allometric scaling enhances stability in complex food webs[END_REF]. This is done because otherwise the number of basal species decreases with increasing C. In this way, [START_REF] Kondoh | Does foraging adaptation create the positive complexity-stability relationship in realistic food-web structure?[END_REF] could obtain positive complexity-stability relations when including foraging dynamics into the model. However, in our investigation a fixation of the number of basal species did not lead to such a positive relation (even though the relation was less negative than before). We therefore used the conventional niche model. The second reason why we did not fix the number of basal species is that this would lead to an exceeding decrease of the fraction of basal species when S is increased, and would thus reduce the energy input to the system. We investigated briefly the effect of fixing the number of basal species in the other food web models. In these models, the number of basal species increases with increasing C, and holding this number constant preserves positive complexity-stability relations, but they are less pronounced. Fixing the number of basal species while S is increased, destroys positive complexity-stability relations where they occurred before. All data shown in this paper are therefore obtained by using the unmodified food web models, where the number of basal species is not fixed by hand.

A c c e p t e d m a n u s c r i p t

Population dynamics

The time evolution of the biomass N i of species i is determined by the balance equation

dN i dt = λ i j∈R i g ij ( N)N i - k∈C i g ki ( N )N k -μ i N i -β i N 2 i (3)
where R i is the set of prey and C i the set of predators of species i. λ i denotes the assimilation efficiency, a value that is different for carnivores and herbivores [START_REF] Brown | Toward a metabolic theory of ecology[END_REF]. For simplicity, we disregard this distinction and set λ i = λ = 0.23. We checked that our results do not change qualitatively when choosing other values of λ. μ i is the mortality of species i, and β i is the intraspecific competition coefficient. Since resources are assumed to have constant size, they are described by

dN 0 dt = dN 1 dt = 0 and N 0 = N 1 = Res . (4)
Simulations have been carried out with Res = 7 or 700.

The energy inflow g ij ( N)N i and outflow g ki ( N )N k contain the functional response g ij , which we chose in our simulations in four different ways:

1. Lotka-Volterra g ij ( N ) = g ij (N j ) = a ij N j .
2. Lotka Volterra with constant foraging efforts,

g ij ( N ) = g ij (N j ) = a ij f ij N j ,
with f ij = 1 B i , and B i the number of prey species of predator species i.

Holling type II

g ij ( N) = a ij N j 1 + k∈C i a ik h ik N k .
4. Holling type II with constant foraging efforts,

g ij ( N ) = a ij f ij N j 1 + k∈C i a ik f ik h ik N k .
Including the factor f ij = 1/B i takes foraging efforts into account. Species with more than one prey divide their available search time among their different prey. Additional parameters in the functional responses are the encounter rate a ij , which was set to a ij = 5 for all links and in all simulations, and handling time h ik , which is the time consumer i needs to deal with and digest one individual of prey k. Simulation results were qualitatively robust against variations of h ik within the tested parameter range between 0.1 and 0.5, so we set the handling time arbitrarily to h ik = 0.3 for all predatorprey interactions. Mortality μ i includes biomass losses due to respiration (metabolism) and other causes, while competition β i limits the growth of a species. This competition may be due to a limitation of nesting sites or territory. We found that a variation of μ i and β i only had a quantitative effect on the food web stability results, but no qualitative effect, as

A c c e p t e d m a n u s c r i p t

long as μ < 0.5a ij and 0.002a ij < βμ < 0.02a ij . We therefore kept their values constant through all simulations and fixed them at μ i = μ = 0.05 and β i = β = 0.4. Due to the intraspecific competition, basal species that feed exclusively on the resources are primary producers with logistic growth, in agreement with other authors [START_REF] Williams | Stabilization of chaotic and non-permanent food web dynamics[END_REF].

With our choice of parameters, interaction rates of species on the lower trophic levels, which had biomasses of the order of 1, were in the saturating regime when Holling type II functional responses were used. Interaction rates of species on the third or fourth trophic level, which typically had biomasses of the order 10 -2 to 10 -4 were on the sloped branch of the functional response. When Holling type II functional response was used with constant foraging efforts, more interaction rates were on the sloped branch, because the inclusion of foraging efforts effectively increases the half saturation density in the functional response of generalist predators.

In our simulations, a species was considered extinct when its biomass dropped below the extinction threshold chosen by us (N i < 10 -6 ), and it was then removed from the system. This could entail the extinction of further species, which depended on the removed species. Such extinction avalanches may lead to the disappearance of a large fraction of the food web.

For each data point, we constructed 200 food webs with the desired topology and values of S and C, and we initiated each web with random biomasses chosen from the interval [0, 1]. The populations of these networks evolved according to the dynamical equations ( 3) and (4), and after 2500 time units we evaluated the average robustness R, i.e., the average fraction of species that have survived. Most extinctions occurred during the first 200 time units, and fixed points or limit cycles are usually reached after no more than 1000 time units. The numerical integration of the dynamical equations was performed using the Runge-Kutta-Fehlberg algorithm with an absolute local error tolerance abs = 10 -5 and a relative local error tolerance rel = 10 -8 .

Simulation results

We first evaluated the robustness of random, cascade, niche and nested hierarchy models by using a Lotka-Volterra functional response. In all cases, we obtained negative complexity-stability relations, i.e. robustness decreased with increasing C or S, if the network remained nontrivial. We call a network nontrivial if it contains several trophic levels, evaluated after population dynamics. The trophic level of a species is defined as the length of the shortest path to the external resources. With larger resources, robustness was higher, but the complexity-stability relations stayed negative.

In the following, we report our results obtained with Holling type II functional response for the different food web models, and compare them to those evaluated with Lotka-Volterra functional response. We devote a separate subsection to each food web model and show the main results for the robustness and the network structure in a set of graphs.

The network structure is specified by the fraction of species in the different trophic levels prior to and after population dynamics. The results show that in almost all cases, including the ones with positive complexity-stability relations, the network structure after population dynamics is nontrivial. The only exception is the cascade model, which is reduced to essentially one level after population dynamics when connectance is high. We also evaluated the network robustness as a function of the final connectance instead of the initial connectance. The behavior remained the same, the results are not shown here.

Apparently, it is more difficult to obtain a positive complexity-stability relation when the species number is varied than when the connectance is varied, since all slopes have a larger value in the latter case. We explain this finding by the fact that the fraction of basal species grows linearly with increasing C, while it remains constant or even decreases with increasing S. This follows directly from the rules for constructing random networks, and it can also be seen in the second row of Figure 1. This means that the energy input into the food web per species grows with increasing C, but stays constant with increasing S, leading to a smaller biomass per species in the second case. The last row of Figure 1 shows that the qualitative shape of these curves is not changed after population dynamics. In particular, the network structure remains nontrivial after population dynamics and even consists of three trophic levels for C < 0.2. However, species on levels higher than the third did not survive.

All model versions that produced negative or neutral complexity-stability relations (B-E) also produced small survival probability results at higher trophic levels. These results are not shown in Figure 1, but are summarized in Table 1, where the fraction of surviving species in the first trophic level and in levels > 1 is given, evaluated for S = 50 and C = 0.15. We chose C = 0.15, because this is an average value of the connectance in empirical food webs [START_REF] Dunne | Food web structure and network theory: the role of connectance and size[END_REF]. (Note that the foodweb data used for our evaluation of vulnerability distributions shown in Figure 6 below include only a subset of the foodwebs in [START_REF] Dunne | Food web structure and network theory: the role of connectance and size[END_REF], creating the wrong impression that the number of links per species does not increase with increasing S.)

We ascribe the smaller robustness of Lotka-Volterra simulations to the smaller dynamical stability. Our simulations showed large population oscillations for models with Lotka-Volterra functional response and large resources. When resources are large, the basal species grow fast, and so do their consumers, and then the predators in the third level, etc. The predators on the higher trophic levels then cause a reduction of the population sizes at lower levels, leading to strong oscillations in the population sizes, which drive many species to extinction. These oscillations are prevented in Holling type II models because population sizes saturate when prey is abundant.

The increase in robustness with increasing resource biomass is easily explained by the larger energy input into higher trophic levels, which then can sustain more species and Our results are different from those obtained earlier by [START_REF] Uchida | Relation between complexity and stability in food webs with adaptive behavior[END_REF], where the incorporation of foraging efforts f ij = 1/B i resulted in a reduced robustness, while we find an increased robustness. This is due to the different parameter regime chosen by Uchida and Drossel. With their values of resource size and interaction strengths the effective growth rate of basal species is about one order of magnitude smaller than it is here and (in the case of Holling type II functional response) even the interaction rates between species on the lower trophic levels are far from saturation. The choice of parameter values in [START_REF] Uchida | Relation between complexity and stability in food webs with adaptive behavior[END_REF] is such that including foraging efforts f ij = 1/B i always decreases all interaction rates which finally leads to a negative growth rate of the species on higher levels even in the absence of predators. In the simulations presented here, where biomass flow to higher levels is large enough, including foraging efforts has a stabilizing effect, because predators with many prey species do not have a larger growth rate than predators with few prey. This prevents the dominance of a single or a few species.

We conclude this subsection by connecting our results to those by May cited at the beginning of this paper. First, we want to point out that our study is fundamentally different, since we did not define stability as the local stability of a dynamical fixed point. Instead, we focussed on species survival and network robustness. Nevertheless, it appears that May's stability criterion is relevant even for our robustness study: We found that foraging efforts f ij = 1/B i need to be included in order to get positive complexitystability relations. This is in agreement with May's stability criterion. Since

f ij 1 S•C , we have α = γ • a ij f ij = γ • a ij f ij = γ • a ij 1 C • S . (5)
with γ = 1+λ 2 • 1 β . This scaling factor takes into account that half of the elements in the community matrix are multiplied with λ and that May scaled his equations in such a way that β = 1. With this, May's criterion becomes

1 > α • (C • S) 1/2 = γ • a ij /(C • S) 1/2 , (6)
which is fullfilled for sufficiently large C or S. This means that the incorporation of foraging efforts inverts May's complexity-stability relation even in the absence of foraging dynamics.

The inclusion of foraging efforts is commonly accepted as being more realistic than using unmodified encounter rates. This means that there is no contradiction at all between May's mathematical result and the finding that ecosystems with many links or many species are stable.

Cascade model

A c c e p t e d m a n u s c r i p t

We can understand this increase in robustness by considering the differences in the structure of the cascade and random models. The main difference is that predators in the cascade model may only establish links to species that have a lower species index, which means that the maximum number of possible links is only half as large (S(S -1)/2) as in the random model and that the number of basal species is much higher in the cascade model. Figure 2 shows that the fraction of basal species is indeed almost twice as high as in the random model (Figure 1). With increasing C, the structure of the network simplifies to the extent that almost all species are basal species, which of course enhances robustness, but makes the network trivial. In order to explain the form of curve C, which starts at considerably lower robustness than curve B and finally bends up, we focus on the attack rates g ij /N j . Let AR B and AR C denote the attack rates for versions B and C respectively. They are

AR B = a ij f ij 1 + B i a ij f ij h ij N j = a ij B i • (1 + B i a ij h ij N j • 1 B i ) , ( 7 )
and

AR C = a ik 1 + B k a ik h ik N k . ( 8
)
for the two models. As long as C is small, the sums in the denominators contain only few terms, and AR B is smaller than AR C because of the factor B i in the denominator of (7). Since lower attack rates lead to a smaller energy outflow from prey species, version B has the higher stability. For larger C, the sum in the denominators will most likely contain the resource biomass, which then becomes the dominant contribution in (8).

Since version C has the larger resource biomass, the denominator of the version C attack rate increases and version C becomes as robust as version B.

In addition to the larger proportion of basal species, other features of the cascade model may also contribute to the enhanced robustness compared to the random model. The cascade model has more top predators than the random model. This is similar to the nested-hierarchy model, which also generates many top predators together with a strong positive slope of the robustness for increasing C. Moreover, in contrast to the random model, the cascade model has no loops, which often lead to the extinction of species on high trophic levels that are preyed upon by species on lower levels.

Similarly to the random model, a positive complexity-stability relation for varying S is obtained only for sufficiently large resources and for f ij = 1/B i (Figure 2). We found that the slope of the complexity-stability relation increases slightly with higher connectance values (results not shown).

Niche model

We could not obtain a positive complexity-stability relation for the niche model even when using foraging efforts, Holling type II functional response, and large resources (Figure 3). This must be due to the specific structural properties of niche model food webs. In fact, the niche model already caused problems when trying to obtain positive complexity-stability relations by including adaptive foraging dynamics [START_REF] Brose | Comment on "Foraging adaptation and the relationship between food-web complexity and stability[END_REF][START_REF] Kondoh | Does foraging adaptation create the positive complexity-stability relationship in realistic food-web structure?[END_REF]. Only when the number of basal species was artificially kept constant, could a positive relationship between C and R be found (Kondoh, 2003b;[START_REF] Kondoh | Does foraging adaptation create the positive complexity-stability relationship in realistic food-web structure?[END_REF][START_REF] Uchida | Relation between complexity and stability in food webs with adaptive behavior[END_REF]. Under no circumstances could a positive

A c c e p t e d m a n u s c r i p t

relationship between S and R be found. When keeping the number of basal species constant or even when increasing it artificially for increasing C, we still obtained negative complexity-stability relations, even though the robustness was increased compared to the original niche model. We found an important difference between the niche model and the cascade and random models when we investigated the distribution of the number of predator species (the so-called vulnerability distribution) for different values of the complexity parameters C and S. Figure 4 shows these distributions for the niche model. The vulnerability distributions of the niche model change qualitatively when increasing connectance. They become broader and develop a hump at large numbers of predator species. This means that with increasing complexity, a larger proportion of species has many predators. The same effect occurred when S was increased (not shown). The vulnerability distributions after population dynamics show that the species with many predators do not survive (Figure 5).

In order to test whether the vulnerability distribution of the niche model is responsible for the negative complexity-stability relation, we ran simulations with a modified niche model: a network was created according to the rules of the niche model, then the links of each species to its prey were rewired to other, randomly chosen species. In this way, the number of prey species of each predator remained the same, but the distribution of the number of predator species became like that of a random model. Simulations with this modified model gave a positive complexity-stability relation under the same conditions as for the random model (i.e., with foraging efforts, with large resources, and with Holling type II functional response). In contrast, a randomization of the distribution of the number of prey species had no effect on complexity-stability relations.

Similarly, a randomization of the connections that preserves the number of predators (and optionally that of prey) of each species, could not change the sign of the complexitystability relation. This last modification removes the unrealistic feature of the niche model that all prey of a species lay within a closed range of niche values. We conclude that while the randomized vulnerability distribution is not necessarily more realistic than the original bimodal one, our results clearly demonstrate that the latter decreases food web robustness when S or C is is increased. These findings beg the question whether the vulnerability distributions generated by the niche model are realistic. In order to answer this question, we evaluated the vulnerability distributions of real food webs. The vulnerability distribution is usually not analyzed and compared to empirical data when the predictive quality of simple topological food web models like the niche model is assessed. [START_REF] Williams | Simple rules yield complex food webs[END_REF][START_REF] Williams | Success and its limits among structural models of complex food webs[END_REF] had compared several niche model features with other models and with real food webs, concluding that the niche model agrees best with data from real webs. (However, these authors had not included models generated by a dynamical process, such as [START_REF] Rossberg | An explanatory model for food-web structure and evolution[END_REF] and [START_REF] Rossberg | Food webs: Experts consuming families of experts[END_REF].) Among the features evaluated by them is the standard deviation of the vulnerability V i of species i, which is the normalized predator count

V i = 1 L/S j a ij . ( 9
)
This quantity does not contain information about the shape of the vulnerability distribution. The bimodal shape has not been seen either in an analytical evaluation

A c c e p t e d m a n u s c r i p t

of the vulnerability distribution in the niche model [START_REF] Stouffer | Quantitative patterns in the structure of model and empirical food webs[END_REF], which was performed in the limit of large species numbers (S ∼ 1000) and a small connectivity z = 2 • L/S = 2 • S • C, (with z 5). For such extreme parameter values (which correspond to a connectance C = 0.005), exponentially decreasing distributions of the numbers of predator and prey species are obtained, which agree well with empirical distributions (which, however typically have connectance values around 0.15, as mentioned before [START_REF] Dunne | Food web structure and network theory: the role of connectance and size[END_REF]). We compare in the following the vulnerability distributions obtained from empirical food web data with those of niche model webs with identical complexity parameters S and C. We analyzed the data of 16 food webs from several terrestrial, freshwater and marine communities (Table 2). The examined food webs have a smaller connectance when S is large, and the product Z = C • S does not show any trend with increasing food web size (Figure 6). In order to improve statistics when evaluating the vulnerability distribution, we divided these food webs into only two classes, "small" and "large", according to the number of species. A food web was considered to be large when S > 50, and small otherwise. Next, we evaluated vulnerability distributions of each food web, accumulated the data and finally normalized the resulting distributions for small and large networks (Figure 7). The results show that an increase in species number leads to a more strongly pronounced exponential distribution of predators per species, a fact that is not mirrored in the niche model. Neither the exponential shape nor the narrowing with increasing S (not even for constant C • S, results not shown) is seen in the niche model vulnerability distribution. We conclude that the peculiar shape of vulnerability distribution of the niche model is an undesired feature of this otherwise realistic food web model. A striking feature of Figure 3 (top curves) is that the curves are strongly clustered according to whether foraging efforts were accounted for (f ij = 1/B i ) or disregarded (f ij = 1). The robustness for the three cases with foraging efforts (versions A,B and E) is by a factor 2 to 3 larger than for the remaining two versions C and D. Some degree of clustering is also visible in the curves R vs C of the cascade model (Figure 2) when C is small, i.e., when the network is nontrivial. We ascribe the strong clustering in the niche model to the fact that the fraction of basal species shrinks with increasing C, i.e. the total number of basal species decreases, while the other models have more links to external resources when C is larger. Additionally, niche model food webs are the only ones in which the proportion of species on level 3 does not decrease when C or S increases. This means that with increasing complexity, the proportion of species on higher levels increases in the niche model. Since species on the second and third level do not feed on the resources, the denominator in the attack rate Eq. ( 8) cannot become large, and a limitation of the attack rates can only be achieved by including foraging efforts, thus increasing the robustness. In contrast, the resource size or the choice of the functional response have no significant effect on the robustness of niche model food webs. random models than those of the niche model. But in contrast to the cascade and random models, the robustness never increased with increasing S. This difference may not be important, as the slope for version A is close to zero in all three models. Just as for the cascade model, version B and C show a (small) positive slope of R vs C. Inspired by our findings for the niche model, we evaluated the distributions of predator and prey numbers in the nested hierarchy model. Figure 9 shows the distribution of predators and prey in the nested hierarchy model for several values of C and S. An increased C leads to less species without predators and to less species with only one prey, whereas a larger value of S leads to more species without predators and to more species with only one prey. A comparison with the experimental data (Figure 7) confirms that larger food webs have more top predators.

Nested Hierarchy model

When randomizing the vulnerability distribution, we obtained positive relations for R as function of S, just as in the niche model. However, since the vulnerability distribution generated by the nested-hierarchy model appears as realistic as other features of this model, the randomization procedure makes the model less realistic and therefore has to be rejected. Additional information about top predator distributions, including their biomasses, in natural food webs would be needed in order to verify to some detail whether the vulnerability distribution of nested hierarchy model is realistic. As top predators represent a considerable threat to prey species, qualitative differences in the predator abundances could have a strong impact on prey population dynamics and complexitystability relations. Based on our present knowledge, we have to conclude that the most realistic food web model has no positive complexity-stability relation when R is evaluated as function of S.

Conclusion

We investigated the robustness and network structure of the four most common food web models (random, cascade, niche and nested hierarchy) in order to establish the conditions under which positive complexity-stability relations emerge, thus allowing for the existence of complex stable foodwebs.

We have shown that for random models large resources together with a functional response of type Holling II and foraging efforts f ij = 1/B i are necessary in order to obtain positive slopes for the robustness as function of C or S. In the cascade model, (weakly) positive relations R vs C are already obtained when the first or third condition is relaxed, while all three conditions are required for a positive slope of R vs S. We argued that our results for random model are consistent with the May criterion and that the larger proportion of basal species is responsible for the higher robustness of the cascade model. The niche model did not produce positive complexity-stability relations. We found that the niche model has an unrealistic distribution of predator numbers per prey species, with a hump at large predator numbers that becomes more pronounced with increasing C or S. This leads to a higher extinction probability with increasing food web complexity. In contrast, empirical data used in this paper and recent studies by other authors [START_REF] Stouffer | Quantitative patterns in the structure of model and empirical food webs[END_REF] suggest that realistic food web models should have exponential distributions of the numbers of prey and predator species. When we randomized the vulnerability distribution in the niche model, we found positive complexity-stability relations under the condition that Holling type II functional response and f ij = 1/B i is used in the simulations. We conclude that the niche model has an unintended unrealistic feature, which has a large effect on the population dynamics. The nested hierarchy model has similar complexity-stability relations as the cascade model when C is increased, but the relation is always negative when S is increased.

Since the distributions of the numbers of prey and predator species seem to agree with empirical data, we cannot blame this negative relation onto unrealistic structural features -even though a randomization of vulnerability distributions leads to positive complexitystability relations.

Out of the four models investigated in this paper, the nested hierarchy model appears to be the most realistic one. It raises, however, the question of why complexity-stability relations with varying S are not positive. There are several possible explanations, including that robustness is not the best measure of stability or that adaptive behavior or some other feature of real ecosystems needs to be included. Or, the nested hierarchy model may also posses structural features that are unrealistic, but that have not yet been uncovered. It would be interesting to compare simulation results for the nested hierarchy model to future empirical data concerning top predators in food webs with the same value of C, but with a different number of species.

Our findings can be summarized as follows: Robustness is larger when Holling type II functional response is used compared to Lotka-Volterra functional response, because oscillations are smaller. Apart from this, the robustness of model food webs is mainly determined by energetic constraints. First, and most obvious, the resources have to be large enough to permit persistence of species on the highest trophic levels. Second, the input rate of energy per species into the food web has to be large enough. This rate is constrained by the fraction of basal species. We hypothesize that most of the robustness differences between the four food web topologies investigated can be explained by this. For increasing connectance, robustness is highest in the cascade model, which also has the highest fraction of basal species, and lowest in the niche model, which is the only one with a decreasing fraction of basal species. Third, the energetic resources of species one the same trophic level should not vary too much. When generalist predators on high trophic levels can feed on each of their prey species with the same rate as specialized predators feeding on only one prey, the latter are outcompeted. This is corroborated by the finding that nearly all simulation runs with f ij = 1/B i yielded higher robustness results than those without partitioning of the foraging efforts. This is to some extend consistent with May's criterion for local stability of random food webs but it also demonstrates that May's criterion does not explain everything. In the cases where we found negative complexity-stability relations in spite of decreasing mean interaction strength as, for example, in the niche model, additional mechanisms like unfavorable topological characteristics have to be invoked to explain the low robustness.

We also investigated the effect of allometric scaling [START_REF] Brose | Allometric scaling enhances stability in complex food webs[END_REF][START_REF] Yodzis | Body Size and Consumer-Resource Dynamics[END_REF] on the network robustness. In agreement with [START_REF] Brose | Allometric scaling enhances stability in complex food webs[END_REF][START_REF] Yodzis | Body Size and Consumer-Resource Dynamics[END_REF], we found that the inclusion of body-size effects generally increases robustness, however we did not find that it can invert the sign of the complexity-stability relation.

In the present study, we analyzed the population dynamics in food webs with constant interaction strengths. An interesting complementing analysis might be to include realistic distributions of link strengths, which are known to include many weak links. (This permits to have many links per species, while at the same time May's criterion can be fulfilled due to a small average coupling strength.) It also remains to be seen how our Table 1: fraction of species surviving on trophic level 1 (TL 1) and higher levels (TL > 1) for versions A to E (defined as in Figure 1).

  The simulation results obtained for random networks are shown in Figure1. A clearly positive complexity-stability relation was only obtained when using Holling type II functional response, large resources (Res = 700), and foraging efforts f ij = 1/B i (version A). When one of these three conditions was relaxed, the relation became neutral (versions B,C for varying connectance, left diagram) or negative (right diagram, for varying species number). With small resources (Res = 7) and without foraging efforts (which is equivalent to setting f ij = 1), the robustness versus connectance relation became negative (version D). As mentioned before, the Lotka-Volterra simulation (version E) gave a negative relation. Similarly, when Lotka-Volterra functional response was used either in combination with small resources or with no foraging efforts, the slope was negative (not shown).

Figure 2

 2 Figure2shows our results for the robustness and the food web structure for the cascade model. The trends of the robustness curves are in general similar to those obtained for random networks. The robustness is for most parameter values higher than in random models, and the curves for versions B (Holling type II with foraging efforts and Res = 7) and E (Lotka-Volterra) have now a slightly positive slope for C > 0.15. A striking feature of Figure2is the behavior of the curve for version C (Holling type II with large resources and without foraging efforts), which bends upwards for C > 0.25.

Finally

  , we evaluated the complexity-stability relations for the nested hierarchy model. The results are shown in Figure8. The data resemble more those of the cascade and

  by including foraging dynamics. Since foraging dynamics is known to increase robustness(Kondoh, 2003a), it may lead to positive complexity-stability relations in some cases where we found negative relations without adaptive foraging.

Figure 1 :

 1 Figures /Captions

Figure 2 :

 2 Figure 2: Cascade model: complexity-stability relations and species distribution on trophic levels for fixed S = 50 and varied C (left) and fixed C = 0.15 and varied S (right). Top: Robustness R. A: functional response Holling II, f ij = 1/B i and Res = 700, B: Holling II, f ij = 1/B i and Res = 7, C: Holling II, f ij = 1 and Res = 700, D: Holling II, f ij = 1 and Res = 7, E: Lotka-Volterra,f ij = 1/B i and Res = 700. Center : initial species distribution on trophic levels (TL) 1,2,3 and > 3, for version A. Bottom: species distribution on trophic levels after population dynamics for version A.

Figure 3 :

 3 Figure 3: Niche model: complexity-stability relations and species distribution on trophic levels for fixed S = 50 and varied C (left) and vice versa (right), C = 0.25. Top: complexity-stability relations. A: functional response Holling II, f ij = 1/B i and Res = 700, B: Holling II, f ij = 1/B i and Res = 7, C: Holling II, f ij = 1 and Res = 700, D: Holling II, f ij = 1 and Res = 7, E: Lotka-Volterra,f ij = 1/B i and Res = 700. Center / Bottom: see previous figures.

Figure 4 :Figure 5 :

 45 Figure 4: Predator distribution obtained from 1000 niche model webs with S = 50 and C = 0.1 (graph A)/C = 0.25 (graph B). The shape of the distribution changes qualitatively when complexity is increased.

Figure 6 :Figure 7 :

 67 Figure 6: Connectivity Z = C • S for empirical food web data vs number of species S. Dotted line: average connectivity.

Figure 8 :

 8 Figure 8: Nested hierarchy model: complexity-stability relations and species distribution on trophic levels for fixed S = 50 and varied C (left) and vice versa (right), C = 0.15. Top: complexity-stability relations. A: functional response Holling II, f ij = 1/B i and Res = 700, B: Holling II, f ij = 1/B i and Res = 7, C: Holling II, f ij = 1 and Res = 700, D: Holling II, f ij = 1 and Res = 7, E: Lotka-Volterra,f ij = 1/B i and Res = 700. Below / Bottom: see preceding figures.

Figure 9 :

 9 Figure 9: Distributions for prey and predator in nested hierarchy model. top: S = 50 and C = 0.1: prey (prey a) and predator (pred a) distribution (thick lines), C = 0.25: graphs indexed with b (thin lines). Bottom: C = 0.15 and S = 20: graphs indexed with c (thick lines), S = 80: graphs indexed with d (thin lines).

Table 2 :

 2 The empirical food webs used for evaluating numbers of predator species, with their connectance C, and species number S. The double line divides the "small" from the "large" food webs.

	Food web data set	Connectance C Food web size S
	Skipwith Pond	0.315	25
	Benguela Stream	0.240	29
	Coachella Valley	0.312	29
	Chesapeake Bay	0.071	31
	St. Martin Island	0.116	42
	St. Marks Seagrass	0.096	48
	Reef	0.220	50
	Grassland	0.026	61
	Bridge Brook Lake	0.22	79
	Ythan Estuary 1991 0.057	83
	Scotch Broom	0.031	85
	Little Rock Lake	0.118	92
	Canton Creek	0.067	102
	Stony Stream	0.070	109
	Ythan Estuary 1996 0.038	124
	El Verde Rainforrest 0.063	155

Acknowledgement

We thank U. Brose for helpful suggestions and for providing the empirical data.


Bogdan Kazmierczak

Tomasz Lipniacki

Regulation of kinase activity by diffusion and feedback

Keywords: reaction-diffusion system, signal transduction, positive feedback, kinase activation, membrane receptors, protein motilitity

In living cells proteins motilities regulate the spatiotemporal dynamics of molecular pathways.

We consider here a reaction-diffusion model of mutual kinase-receptor activation showing that the strength of positive feedback is controlled by the kinase diffusion coefficient. For high diffusion, the activated kinase molecules quickly leave the vicinity of the cell membrane and cannot efficiently activate the receptors. As a result, in a broad range of parameters, the cell can be activated only if the kinase diffusion coefficient is sufficiently small. Our simple model shows that change in the motility of substrates may dramatically influence the cell responses.

Introduction

Regulatory networks process cellular signals in time and space enabling cell self-organization (see [START_REF] Kholodenko | Cell-signalling dynamics in time and space[END_REF][START_REF] Karsenti | Self-organization in cell biology: a brief history[END_REF] for reviews). The temporal dynamics is coupled with spatial gradients of concentrations or activity. For example, kinase cascades can emerge from receptors and transmit signals from the cell membrane to the nucleus. In this case the gradient of active kinase activity develops since phosphorylation and dephosphorylation proceed at different cellular locations, respectively cell membrane and cell volume. Due to the estimations of [START_REF] Brown | Spatial gradients of cellular phospho-proteins[END_REF], basing on measured values of protein diffusion coefficients and phosphatase activities, gradients of kinase activity are potentially very large. In a simple system in which kinase molecules are phosphorylated at the cell membrane and dephosphorylated by a phosphatase molecules located homogeneously in the cell cytosol (analyzed by [START_REF] Brown | Spatial gradients of cellular phospho-proteins[END_REF] small diffusion implies high gradient and low kinase activity in the cell center. The problem of receptor-kinase interaction has been also studied in the context of diffusion with obstacles in the stochastic numerical simulations of bacterial chemotaxis [START_REF] Lipkow | Simulated diffusion of phosphorylated CheY through the cytoplasm of Escherichia coli[END_REF]. One of the conclusions of [START_REF] Lipkow | Simulated diffusion of phosphorylated CheY through the cytoplasm of Escherichia coli[END_REF] is that crowding results in a fall of the apparent diffusion coefficient and at the anterior end, where CheY is phosphorylated, the local concentration of CheYp increases and therefore accelerates the response of the anterior close motor. At the other, posterior, end of the cell, the local CheYp concentration is reduced by the need to diffuse through the obstacles and the responses of motors in this region is consequently delayed.

A c c e p t e d m a n u s c r i p t

Here we consider a similar model to the one analyzed by [START_REF] Brown | Spatial gradients of cellular phospho-proteins[END_REF], but assume the mutual receptor -kinase activation. Membrane receptors can bind extracellular ligands, that leads to cascade of molecular processes inside the cell and formation of the active receptor complex. In many cases, receptor activation requires phosphorylation. Almost all G-protein coupled receptors (GPCRs) are regulated by phosphorylation, see [START_REF] Tobin | G-protein-coupled receptor phosphorylation: where, when and by whom[END_REF] for review. Engagement of immunoreceptors (TCR, BCR, FcR) leads to activation of different members of the Src kinase family, which includes Lck (for T-cell, [START_REF] Housden | Investigation of the kinetics and order of tyrosine phosphorylation in the T-cell receptor ξ chain by the protein tyrosine kinase Lck[END_REF], Fyn and Lyn (for B and mast cells, [START_REF] Gauld | Src-family kinases in B-cell development and signaling[END_REF]. Src kinases then phosphorylate immunoreceptor tyrosine-based activation motifs (ITAMs) contained within the immunoreceptors themselves or in receptor-associated molecules, see [START_REF] Abram | Convergence of immunoreceptor and integrin signaling[END_REF] for review. This may lead to positive feedback, in which active receptors send signal to kinase and in turn are activated by the same kinase species or by one of the downstream kinases. In this study we consider the simplest situation, in which receptors are activated by the same kinase species they activate.

We will show, that in the case of mutual receptor -kinase activation, in a broad range of parameters controlling the process, the cell becomes activated only if the kinase diffusion is sufficiently small. For large diffusion, the activated kinase molecules quickly leave vicinity of cell membrane, and the positive feedback coupling kinases with receptors becomes inefficient.

Model formulation

As said, we will assume that membrane receptors bind extracellular ligands, that leads to a cascade of processes and receptor activation. At constant extracellular cytokine concentration, a steady state uniform surface concentration of ligand-bound receptors P = const is established. We will consider two cases. Firstly, following [START_REF] Brown | Spatial gradients of cellular phospho-proteins[END_REF] we will assume that all ligandbound receptors are active. Secondly, we will assume that the limiting step in the formation of the active receptor complex is its phosphorylation by the kinase. In turn, active receptors may activate kinase molecules, that defines the positive feedback in the regulation process. The activated kinase may freely diffuse over entire cell volume, where they are inactivated by uniformly distributed phosphatases.

The cell will be modeled geometrically as a ball B(0, r 0 ) of radius r 0 , centered at the origin of the coordinate system. We restrict to the spherically symmetric case and we will use the following notation : The active kinase concentration satisfies

∂K ∂t = d 1 ∇ 2 K -b 1 K, (1)
where b 1 > 0 is the kinase dephosphorylation rate due to the action of uniformly distributed phosphatases. The flux of the active kinase results from its phosphorylation by the surface receptors implying the Robin type boundary condition,

Φ K = a 1 R (Q -K b ) = d 1 n • (∇K) b , (2)
where n is a unit vector normal to cell surface and subscript b denotes the boundary value for (r = r 0 ).

As already mentioned, regarding the receptors activation we will consider two different cases;

(1) Activity of the receptors is independent of intracellular processes, but is controlled by binding and dissociation of some extracellular ligand (present at the constant concentration). In such a case we may assume that all the bound receptors are persistently active i.e. R(t) = const = P .

(2) Alternatively, we may assume that the limiting step in the receptor activation is its phosphorylation by the kinase, that defines the positive feedback in receptor-kinase activation,

dR dt = a 2 K b (P -R) -b 2 R. (3)
In the further consideration we will assume that all the reaction rate coefficients, a 1 , a 2 , b 1 , b 2 , and diffusion constant c 1 are positive. In the non-dimensional units τ = t b 1 , ρ = r/r 0 , Eq. 1 reads

∂K ∂τ = d∇ 2 K * -K * , (4)
where d = d 1 /(b 1 r 2 0 ) plays the role of the non-dimensional diffusion coefficient and K * = K/Q. We may thus rewrite Eqs. 2 and 3 as

aR * (1 -K * b ) = dn • (∇K * ) b , (5)
A c c e p t e d m a n u s c r i p t

dR * dτ = qK * b (P * -R * ) -bR * , (6)
where

a = a 1 P/b 1 , R * = R/(r 0 Q), q = a 2 Q/b 1 , b = b 2 /b 1 and P * = P/(r 0 Q).
From now on, for the sake of simplicity all the asterisks will be omitted. Let us notice that ρ ∈

[0, 1], K ∈ [0, 1],
R ∈ [0, P ] and a, b, d, q and P are real and nonnegative.

Results

Limit of infinite diffusion

d → ∞, K = K(τ ), R = R(τ ).
For the infinite diffusion, the active kinase concentration is uniform, and the system of Eqs. 4-6 is equivalent to the system of two ordinary equation, for K(τ ), R(τ ),

dK dτ = 3aR(τ ) (1 -K(τ )) -K(τ ), (7)
dR dτ = qK(P -R(τ )) -bR(τ ). (8)
Formally, we can obtain the above system by integrating Eq. 4 over the ball B(0, 1) and using the Gauss theorem. Let us note that the compact region D := [0, 1] × [0, P ] is invariant with respect to the flow generated by the above system. That is to say, if {K(0), R(0)} ∈ D, then for arbitrary

τ ≥ 0, {K(τ ), R(τ )} ∈ D. System 7-8 has two steady states: {K 1 , R 1 } = {0, 0} and {K 2 , R 2 },
where

K 2 = 3qaP -b q(3aP + 1) , R 2 = 3qaP -b 3a(b + q) . ( 9
) For b < 3qaP the steady state point {K 2 , R 2 } is stable, while the point {K 1 , R 1 } is unstable, for b ≥ 3qaP , K 2 < 0, R 2 < 0 and the point {K 2 , R 2 } is unstable, while the {K 1 , R 1 } is stable.
In other words, restricting to subdomain D, the system has one stable steady state {0, 0}, for b ≥ 3qaP , or {K 2 , R 2 } for b < 3qaP .

In the case without the feedback, (i.e. when R(τ ) = P ) Eq. 7 can be solved analytically,

A c c e p t e d

m a n u s c r i p t

K(τ ) = K(0) - 3aP 3aP + 1 exp[(1 + 3aP )τ ] + 3aP 3aP + 1 (10)
and has unique stable state K 3 = 3aP/(1 + 3aP ).

Finite diffusion -steady state analysis

(1) Case without feedback

In spherical coordinates Eq. 4 reads

∂K ∂τ = d 1 2 ∂ ∂ 2 ∂K ∂ , (11)
and has unique steady state solution

K( ) = K c (e α -e -α ) 2 α , (12)
where α = d -1/2 and K c = K(0). Using Eq. 5 (with R = P ) and Eq. 12 we may calculate boundary value

K b = K(1) K b = aP α 2 (e 2α -1) 1 + α -aα 2 P + e 2α (α + aα 2 P -1)) , (13)
and then

K c = 2K b α/ e α -e -α = 2aP α 3 e α 1 + α -aα 2 P + e 2α (α + aα 2 P -1)) . ( 14
)
In the limit of infinite diffusion α → 0, the active kinase distribution is uniform, K( ) ≡ K c , with K c = K 3 . The total amount of active kinase K tot is

K tot = 4π 1 0 2 K( )d = 4πK c αcosh(α) -sinh(α) α 3 . ( 15
)
Differentiating K b , K c and K tot several times with respect to α, one can show that for a > 0 and P > 0:

1) K b is a monotonically increasing function of α = d -1/2 ,
2) K c and K tot are monotonically decreasing functions of α.

A c c e p t e d m a n u s c r i p t

Large diffusion enhances the flux of the active kinase from the cell membrane towards the cell center. Thus, as one could expect, both K c and K tot are growing functions of diffusion coefficient.

Simultaneously, the active kinase concentration close to the boundary decreases with increasing diffusion, Fig. 1.

(2) Case with feedback

In the case with feedback the spatial profile of K( ) is the same as in the case without feedback;

Eqs. 5-6 gives us two values of K b (or

K c = K b * 2α/ (e α -e -α
)) corresponding to one stable and one unstable steady state solution

K c1 = 0, K c2 = 2αe α be 2α (1 -α) -b(α + 1 + qaP α 2 e 2α -1 ) q (e 2α -1) (1 + α -aα 2 P + e 2α (α + aα 2 P -1)) . ( 16
)
The stable solution K( ) is given by

K( ) = K c (e α -e -α ) 2 α . ( 17
)
where K c = max(K c1 , K c2 ), i.e. for K c2 > 0 the stable solution is positive, while for K c2 < 0, K( ) ≡ 0. The global stability (with respect to a perturbation without spherical symmetry) of K( ) solution is proved in the Appendix.

Let us note, that in the limit of the infinite diffusion coefficient d → ∞, i.e. α → 0, K( ) ≡ K c = K 2 and thus one obtains the same solution as given in Eq. 9. In the opposite limit d → 0 (i.e. α → ∞), K b = 1, K c = 0, K tot = 0. In further analysis we set q = a = P = 1, and consider the steady state kinase activity profiles with respect to two nondimensional parameters: α = d -1/2 , and b (receptor inactivation constant).

In Fig. 2 we analyze the dependence of the stable steady state K( ) on the diffusion parameter For small b, K c decreases with growing α (as in the case without feedback), however for larger b, K c (α) has a maximum for some α m (b) > 0. The existence of such an "optimal" α m is due to interplay of two counter-effects:

α = d -1/2 . With
(1) large diffusion (small α) speeds translocation of active kinase, so they have a larger chance to remain phosphorylated until they reach the cell center,

(2) simultaneously large diffusion attenuates the positive feedback coupling kinases with receptors.

In contrast to the case without feedback, the total amount of active kinase K tot

Discussion

Dynamics of molecular pathways is determined by both, chemical reaction rules and localization of substrates that in turn is governed by diffusion or transport. We considered here a simple theoretical model of mutual receptor-kinase activation in which the kinase molecules, are phosphorylated by the receptors at the cell membrane and may freely diffuse in the cell volume, where they are dephosphorylated with time-and space-independent dephosphorylation rate. The positive feedback, considered in the model, arises since activated kinase may in turn activate receptors.

In the case without feedback, active kinase concentration in the cell center and the total amount of active kinase are monotonically growing functions of diffusion coefficient. However, the presence of the positive feedback causes that the concentration of the active kinase in the cell volume is a nontrivial function of the diffusion coefficient. The active kinase concentration in the cell center depends on the two opposite, diffusion controlled effects: the kinase activity profile along the cell

A c c e p t e d m a n u s c r i p t

radius is flatter for large diffusion, and the strength of positive feedback controlling active kinase concentration at the boundary is stronger for small diffusion. As a result in a broad range of parameters a, q, P and b controlling mutual kinase-receptors activation and inactivation rates the maximum value of the active kinase concentration in the cell center is reached for some finite value of the kinase diffusion coefficient. Moreover, for the large receptor inactivation rate (b ≥ 3qaP ) the active kinase concentration is everywhere zero for the infinite diffusion, but it is positive for the sufficiently small diffusion. Interestingly, for constant a, q and P , even for arbitrarily large inactivation coefficient b, there exists a positive steady state solution for the active kinase concentration if the diffusion is sufficiently small.

In living cells the diffusion and thus spatiotemporal localization of substrates can be controlled in a number of ways; Molecules can bind to a larger molecules of lower motility called buffers, or to the cell membrane and other structural elements directly or with help of the, so called, anchoring proteins. On the cell membrane receptors can form larger complexes of lower motility, or get localized within lipid rafts. Relevant to our model, cell membrane can create microdomains which trap signalling molecules, like Lck kinase, that activates TCR receptors, [START_REF] Douglass | Single-molecule microscopy reveals plasma membrane microdomains created by protein-protein networks that exclude or trap signaling molecules in T cells[END_REF].

Major simplification of our study is that it does not account for macromolecular crowding and presence of diffusion obstacles within the cell (organelle, cellular structures). The macromolecular crowding may have non-trivial effect on molecular association in the cell, possibly increasing its rate by limiting the volume in which molecules are free to diffuse (see [START_REF] Minton | The influence of macromolecular crowding and macromolecular confinement on biochemical reactions in physiological media[END_REF][START_REF] Zimmerman | Macromolecular crowding: Biochemical, biophysical, and physiological consequences[END_REF][START_REF] Zimmerman | Macromolecular crowding: Biochemical, biophysical, and physiological consequences[END_REF].

Our simple model provides an example in which diffusion controls the strength of the feedback regulation and thus the dynamics of kinase activation. In the considered model, for a broad range of parameters the cell can be activated only when the kinase diffusion coefficient is sufficiently small, i.e. when reacting kinases and receptors are well collocalized that enables their mutual activation.

Acknowledgments

This study was supported by Polish Committee for Scientific Research Grants No. 4 T07A 001 30

and No. 1P03A01230.

Figure captions

K tot (α, b) = 0; Δ tot = K tot (α, b) -K tot (0, b) > 0; Δ tot ≤ 0. In domain Δ tot > 0, isolines of Δ tot are
shown. For both plots the remaining parameters are fixed: q = a = P = 1.

Appendix. Stability of spherically symmetric solutions

Let K(ρ) and R = const denote the spherically symmetric stationary solution of system 4-6. For the initial data not coinciding with ( K(ρ), R), the solution (K(x, τ ), R(τ )) to system 4-6 will be in general different, i.e.

K(x, τ ) = K(ρ) + δK(x, τ ), R(x, τ ) = R + δR(x, τ ).
Our task here is to examine the asymptotic in time behavior of the functions K and R. We will confine here to initial data preserving the positivity of the functions K and R. The equations for δR and δK read: Then, the above equations can be written as:

∂δk ∂τ = d∇ 2 δk -(1 -ε)δk (19) dδr dτ = q(P -R)δk b -(q K b + b + qδk b exp (-ετ ) -ε)δr (20) aδr(x, τ )(1 -K b ) = dn • ∇δk + a R + δr(x, τ ) exp (-ετ ) δk b . ( 21
)
The idea of the stability proof is to construct a time independent sub-and supersolution pairs:

(-δr -, δk -) and (δr + , δk + ). As ε > 0 this will prove that δR(x, τ ) and δK(x, τ ) tend to 0 as

τ → ∞. Let δK ± (ρ) = ±δK ±b φ(ρ; αε), (22)
where αε = 1ε d

Let us recall that φ is an increasing function of ρ and φ(1; α) = 1 for all α ∈ (0, ∞). In fact

K(ρ) = K b φ(ρ; α 0 ). (24)
Let

δK -b = K b -ψ(ε) (25)
where ψ(ε) 0 as ε → 0. This function will be specified later. It follows from 23, 24 and the continuity of the function φ(ρ; αε) with respect to the parameter ε that if K(x, 0) > η > 0 for

|x| ≤ 1, we can find ε > 0 so small that K(x, 0) > K(ρ(x)) -δK -b φ(ρ(x); αε) > 0 for all x inside the sphere. Let δK + be at least so large that K(x, 0) < K(ρ(x)) + δK +b (ρ(x)), implying that δK(x, 0) < δK +b (ρ(x)). Let

W R = q(P -R) q K b -qδK -b + b -ε = q(P -R) qψ(ε) + b -ε . (26)
We will assume that

δR ± = δK ±b (W R + ε). (27)
Let us note that as

(P -R)/ R = b/(q K b ), (28)
it follows from Eq.27 that:

δR -= R -C 2 ψ(ε) + C 2 ε + o(ψ(ε)) + o(ε). (29)
for some positive constants C 2 and C 2 . If, for small ε ≥ 0, ψ(ε)/ε ≥ C 3 > 0 with C 3 sufficiently large, then δR -< R. On the other hand, for any R(0) > 0 we can find ε > 0 sufficiently small such that R(0) > R -δR -. Obviously, we may also choose δK + so large that R(0) < R + δR + .

Let

δr ± = δR ± , δk ±b = δK ±b . (30)
Let δk ± (•) be constant in time spherically symmetric solutions to Eq.19:

δk ± (ρ) = ±δk ±b φ(ρ; αε). (31)
We will prove that (-δr -, δk -) and (δr + , δk + ) defined in 30 and 31 are respectively the sub-and supersolution pairs for system 19-21. So, after putting δr = -δr -in Eq.20 we infer that the right hand side is positive if δk b ≥ -δk -b . In the similar way, after putting δr = δr + in Eq.20 we infer that the right hand side is negative if δk b ≤ δk +b . It follows that if δk + ≥ δk(x, τ ) ≥ -δk -at the boundary, then δr(x, τ ) ∈ (-δr -, δr + ). We have thus to prove the corresponding properties of the functions δk -and δk + with respect to Eqs.19-21. First, as we said, the functions δk ± (ρ) satisfy Eq.19. So, to prove that δk -, defined in ( 22) is a subsolution it suffices to show, as (δk - putting δr(x, τ ) = -δr -, we conclude that inequality 32 is implied by the inequality:

(1

-K b + δK -b ) q(P -R) q K b -qδK -b + b -ε + ε < R + a -1 w 1 -ε d
where we used 27. Now, due to 28, (P -R)/ R = b/(q K b ). Hence, we obtain the condition:

b(1 -K b + δK -b )/(q K b -qδK -b + b -ε) + O 1 (ε) < K b + K b R a -1 w 1 -ε d .
As the spherically symmetric solution satisfies the equality

a R 1 -K b = w 1 d K b (33) then b(1 -K b + δK -b )/(q K b -qδK -b + b -ε) + O 1 (ε) < K b + (1 -K b )w 1-ε d w 1 d -1 ≤ K b + (1 -K b )(1 -ν(ε))
for some given smooth function ν. Using 25 we arrive at the inequality b(1ψ(ε))

(qψ(ε) + b -ε) + O 1 (ε) < 1 -(1 -K b )ν(ε).
This inequality can be satisfied if only ψ(ε) ≥ C(ε + ν(ε)) with C > 0 sufficiently large.

In the last step we have to prove that δr(x, τ )(1 -K b ) ≤ R + δr(x, τ ) exp(-ετ ) δk +b + a -1 w 1-ε d δk +b for all δr(x, τ ) ∈ (-δR -, δR + ). As δk +b = δK +b > 0 and δR -< R then δK +b R + δr(x, τ ) > 0 for δr(x, τ ) ∈ (-δR -, δR + ). Thus we have only to prove that δr(x, τ )(1 -K b ) ≤ a -1 w 1-ε d δK +b . Obviously, it suffices to show it for δr = δr + . Let us recall that δr + = δK +b (W R + ε) = δK +b q(P -R) qψ(ε) + bε + δK +b ε.

Hence we have to prove that 4. Figure

q(P -R) qψ(ε) + b -ε + ε (1 -K b ) < a -1 w 1 -ε d .

  r) -the concentration of the active kinase Q = const -the total concentration of the kinase R(t) -the surface concentration of the active receptors P = const -the total surface concentration of the ligand bound receptors (active and inactive) Φ K (t) -the flux of the active kinase

  respect to the receptor dephosphorylation coefficient b, we may distinguish three cases, shown in Panels A, B and C. For small b (Panel A), the dependence of active kinase concentration profiles K( ) on α is similar as in the case without feedback; the larger is α, the steeper is the active kinase profile, with higher value at the boundary and the lower value in the cell center. Qualitatively different is the case shown in Panel B for larger values of b. For b = 2.5, the active kinase concentration is higher across the whole cell for some finite diffusion (α = 2) than for the infinite diffusion (α = 0). In the case of large dephosphorylation parameter (Panel C, b = 4) for infinite, or large diffusion, K( ) ≡ 0, while for smaller diffusion K( ) > 0. This somehow surprising effect is due to the fact that the strength of positive feedback is controlled by the diffusion. For kinase remain longer in the vicinity of the membrane and may activate the receptors more effectively.In Figs.3A and 3Bwe analyze K b (α) and K c (α) for four different values of b. As can be expected , K b is a growing function of α.

  (α) for large b has maximum for some α m (b) > α m (b) > 0. Let us note also, that both α m (b) and α m (b) are growing functions of b, diverging logarithmically to infinity with b. As shown in Fig. 4 there are unbounded parameter domains D c and D tot in (α, b) plane for which respectively Δ c := K c (α, b) -K c (0, b) > 0 and Δ tot := K tot (α, b) -K tot (0, b) > 0. Since for finite diffusion K b (α) > K c (α): D c ⊂ D tot . In the case of the infinite diffusion (α = 0), the positive solutions K(ρ, α, b) > 0 are restricted to domain b < 3. In the case of finite diffusion, for arbitrarily large b there exists such α(b) that K(ρ, α, b) > 0.

Fig. 1 .Fig. 2 .

 12 Fig. 1. The case without feedback. Profiles of active kinase concentration K(ρ) for different values of α = d -1/2 . For all plots the remaining parameters are fixed: a = P = 1.

  ∂δK ∂τ= d∇ 2 δK -δK inside the sphere anddδR dτ = qδK b (P -R) -(q K b + b)δR -qδK b δR aδR(1 -K b )a RδK b -aδRδK b = cn • ∇δK b .Here we used the fact that R and K satisfy system 4-6. For ε > 0 an arbitrarily small positive number, let δk = δK exp ετ , δr = δR exp ετ (18)

  e -α ) (e αe -α ).

  ) |B = w( 1-ε d )(-δk -b ), that for all δr(x, τ ) ∈ (-δr -, δr + ) we haveδr(x, τ )(1 -K b + δk -b ) ≥ -Rδk -ba -1 w 1ε d δk -b .(32)(See, chapter 2 of Pao 1992. Note that the coefficient by δk b in 21 is positive.) Here = φ ,ρ (ρ = 1, s) is smooth and monotonically increasing function of s for s ≥ (2d) -1 . As δk -b = δK -b < K b then,

  Fig. 1

A c c e p t e d m a n u s c r i p t A c c e p t e d m a n u s c r i p t

Using 28 and 33 we obtain as before that we have to satisfy the condition

This is implied by the inequality:

where ν(ε) → 0 as ε → 0. As before, this condition can be satisfied, if only ψ(ε) ≥ C(ε + ν(ε)) with C > 0 sufficiently large. Finally, taking into account what we have shown and using Theorem 2.1.2 from Pao 1992, we come to a conclusion that -δr -< δr(x, τ ) < δr + and δk -(ρ(x)) ≤ δk(x, τ ) ≤ δk + (ρ(x)) for all τ ∈ (0, ∞).

So, due to the definition 18 we conclude that δK(x, τ ) and δR(x, τ ) tend to zero in the supremum norm as t → ∞.

We have thus shown that, if the initial data K(x, 0) and R(0) are positive, then the solution (K(x, τ ), R(τ )) to system 4-6 tends to the unique spherically symmetric solution ( K, R) as τ → ∞.