Jacques Droulez

David Colliaux

Audrey Houillon

Pierre Bessière

Toward Biochemical Probabilistic Computation

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Toward Biochemical Probabilistic Computation

Jacques Droulez a,b, * , David Colliaux a , Audrey Houillon b,1 , Pierre Bessière a,b a CNRS/Sorbonne Universités/UPMC/ISIR, Paris, France b CNRS/Collège de France/LPPA, Paris, France * E-mail: Jacques.Droulez@College-de-France.fr [START_REF] Jaynes | Probability Theory: the Logic of Science[END_REF]

Abstract

To account for the ability of living organisms to reason with uncertain and incomplete information, it has been recently proposed that the brain is a probabilistic inference machine, evaluating subjective probabilistic models over cognitively relevant variables. A number of such Bayesian models have been shown to account efficiently for perceptive and behavioral tasks. However, little is known about the way these subjective probabilities are represented and processed in the brain. Several theoretical proposals have been made, from large populations of neurons to specialized cortical microcircuits or individual neurons as potential substrates for such subjective probabilistic inferences. In contrast, we propose in this paper that at a subcellular level, biochemical cascades of cell signaling can perform the necessary probabilistic computations. Specifically, we propose that macromolecular assemblies (receptors, ionic channels, and allosteric enzymes) coupled through several diffusible messengers (G-proteins, cytoplasmic calcium, cyclic nucleotides and other second messengers, membrane potentials, and neurotransmitters) are the biochemical substrates for subjective probability evaluation and updating. On one hand, the messengers' concentrations play the role of parameters encoding probability distributions; on the other hand, allosteric conformational changes compute the probabilistic inferences. The method used to support this thesis is to prove that both subjective cognitive probabilistic models and the descriptive coupled Markov chains used to model these biochemical cascades are performing equivalent computations. On one hand, we demonstrate that Bayesian inference on subjective models is equivalent to the computation of some rational function with nonnegative coefficient (RFNC), and, on the other hand, that biochemical cascades may also be seen as computing RFNCs. This suggests that the ability to perform probabilistic reasoning is a very fundamental characteristic of biological systems, from unicellular organisms to the most complex brains.

Author Summary

Living organisms survive and multiply even though they have uncertain and incomplete information about their environment and imperfect models to predict the consequences of their actions. Bayesian models have been proposed to face this challenge. Indeed, Bayesian inference is a way to do optimal reasoning when only uncertain and incomplete information is available. Various perceptive, sensory-motor, and cognitive functions have been successfully modeled this way. However, the biological mechanisms allowing animals and humans to represent and to compute probability distributions are not known. It has been proposed that neurons and assemblies of neurons could be the appropriate scale to search for clues to probabilistic reasoning. In contrast, in this paper, we propose that interacting populations of macromolecules and diffusible messengers can perform probabilistic computation. This suggests that probabilistic reasoning, based on cellular signaling pathways, is a fundamental skill of living organisms available to the simplest unicellular organisms as well as the most complex brains.

Introduction

The information available to living organisms about their environment is uncertain, not only because biological sensors are imperfect, but more importantly because sensors inevitably provide an incomplete, partial description of the environment. Moreover, timing is a crucial constraint for biological systems. During a fight-or-flight dilemma, animals must quickly decide to fight or flee, and they can never be sure of the good or bad consequences of their decision. Incompleteness is therefore a key notion for an autonomous agent facing the complexity of the world. There is now growing evidence that probabilistic reasoning is a rigorous and efficient way to represent partial or incomplete knowledge and to answer questions optimally that have no uniquely defined solution [START_REF] Jaynes | Probability Theory: the Logic of Science[END_REF].

Perception is a well-known example of ill-posed problems, because an indefinite number of object characteristics can theoretically give rise to the same set of sensory data. For instance, an indefinite number of objects with various shapes, sizes and movements can induce exactly the same retinal projection. A number of studies have shown that Bayesian models can accurately account for various aspects of perception [START_REF] Knill | Perception as bayesian inference[END_REF][START_REF] Weiss | Motion illusions as optimal percepts[END_REF][START_REF] Rao | Probabilistic models of the brain: perception and neural function[END_REF][START_REF] Mamassian | Bayesian modelling of visual perception[END_REF][START_REF] Ernst | Humans integrate visual and haptic information in a statistically optimal fashion[END_REF][START_REF] Laurens | Bayesian processing of vestibular information[END_REF] and sensory-motor integration [8][START_REF] Körding | Bayesian decision theory in sensorimotor control[END_REF][START_REF] Colas | A unified probabilistic model of the perception of three-dimensional structure from optic flow[END_REF][START_REF] Todorov | General duality between optimal control and estimation[END_REF].

Nonprobabilistic models assume that the brain computes internal estimates of relevant state variables such as motion, object distance, and color. Each variable is supposed to have a unique estimate, on which no evaluation of uncertainty is performed. In contrast, Bayesian models assume that the brain evaluates the probability corresponding to each possible value of the relevant variables. Probability computation results from straightforward application of the Bayes (i.e., multiplicative) and marginalization (i.e., additive) rules, which can been seen as a generalization of logical inference to probability distributions [START_REF] Bessière | Probabilistic Reasoning and Decision Making in Sensory-Motor Systems[END_REF][START_REF] Bessière | Bayesian Programming[END_REF]. Following Pearl [START_REF] Pearl | Probabilistic reasoning in Intelligent Systems: Networks of Plausible Inference[END_REF] and Jaynes [START_REF] Jaynes | Probability Theory: the Logic of Science[END_REF], we define probabilistic reasoning as the ability to perform inference within a probabilistic knowledge base. In the following, we will call a subjective Bayesian model the specification of the variables of interest, their conditional dependencies, the parametric forms, and the way probability distributions can be inferred.

Acknowledging the efficiency of probabilistic reasoning in accounting for a large variety of perceptive reports or motor behaviors, one of the main scientific challenges is to explicitly demonstrate that the brain, and more generally biological systems, can effectively perform probabilistic computation. The problem is to show possible correspondences between subjective Bayesian models and descriptive models of biological systems and signal processing.

Most existing studies proposed assemblies of neuronal cells or single neurons as the appropriate level of analysis to explain how a brain could perform probabilistic inference. Several authors have proposed that the firing rate of a group of cells within a given temporal window [START_REF] Zemel | Probabilistic interpolation of population code[END_REF][START_REF] Deneve | Reading population codes: a neural implementation of ideal observers[END_REF][START_REF] Gold | Banburismus and the brain: decoding the relationship between sensory stimuli, decisions, and reward[END_REF] could represent probability distributions over state variables. The mean firing rate might be well approximated by a graded value, to which a subjective probabilistic meaning can be attributed: for instance, the probability that a given proposal over a state variable is true, or its log-likelihood ratio. Other approaches are based on the Poisson-like variability of spike trains [START_REF] Ma | Spiking networks for bayesian inference and choice[END_REF] or kernel convolution for encoding/decoding spike trains [START_REF] Deneve | Bayesian spiking neurons i: Inference[END_REF].

In contrast, we consider the molecular scale as an adequate framework to solve this matching problem between subjective Bayesian models and descriptive biological ones.

Populations of macromolecules in their various conformational states and diffusible messenger concentrations are assumed to be the substates used at subcellular level to represent and to compute probability distributions. Our proposal is that the biochemical processes involved in cell signaling can perform the elementary computations needed for subjective probabilistic reasoning, and that this biochemical computation is used extensively by the brain as an elementary component: a Bayesian "nanogate." Specifically, we propose that macromolecular assemblies (receptors, ionic channels, and allosteric enzymes) coupled through several diffusible messengers (G-proteins, cytoplasmic calcium, cyclic nucleotides and other second messengers, membrane potentials, and neurotransmitters) are the biochemical substrates for subjective probability evaluation and updating. On the one hand, messenger concentrations play the role of parameters encoding probability distributions, on the other hand, allosteric conformational changes compute the probabilistic inferences. Diffusible messenger concentrations, including electric charge density, control the probability of conformational changes, which are in turn responsible for inflow and outflow rates of messengers and are then controlling their kinetics.

We started exploring these ideas in [START_REF] Houillon | The probabilistic cell: Implementation of a probabilistic inference by the biochemical mechanisms of phototransduction[END_REF], where the kinetics of a rhodopsin channel were shown to converge to a posterior distribution of states in a hidden Markov model, the hidden state representing the presence or absence of light. A similar kinetic scheme derived from Bayesian formulation of receptor activations was compared with other common signaling schemes in [START_REF] Kobayashi | Implementation of dynamic bayesian decision making by intracellular kinetics[END_REF] to emphasize its optimality. Rather than considering a binary (ON or OFF) hidden state, multiple activated states were considered in [START_REF] Siggia | Decisions on the fly in cellular sensory systems[END_REF] (like self and nonself in the immune system) and they discussed how Koshland-Goldbeter could solve this decision problem. Napp and Adams proposed in a recent paper [START_REF] Napp | Message passing inference with chemical reaction networks[END_REF] a procedure to compile a probabilistic graphical model into a chemical reaction network.

Results

Subjective probabilistic models (often called Bayesian models) include variables describing states of the world highly relevant for the organism, such as the presence of food or predators. The values taken by these variables cannot be known with certainty by the organism. However, according to the Bayesian approach, the organism can evaluate the probability distribution over these variables on the basis of specific observations, such as the detection of light or the detection of an odorant molecule.

Subjective probabilistic reasoning is the process by which the probabilities of relevant variables are computed from a set of observations and a set of priors and conditional dependencies between variables [START_REF] Bessière | Bayesian Programming[END_REF].

Descriptive biochemical models, reduce biochemical signaling mechanisms to two sets of interacting substrates: (i) macromolecular assemblies (such as receptors, ionic channels, or activatable enzymes), which have a nearly fixed spatial distribution (in the time scale we considered here), but a variable conformational state distribution; and (ii) diffusible messengers (such as G-protein α-subunits, cyclic nucleotides, cytoplasmic calcium, and membrane potentials) that have a fixed conformational state, but a variable spatial distribution.

The messengers, when they bind to their specific receptor sites on macromolecules, induce conformational changes. The macromolecules, when they are in active states, release or remove messengers, thus changing the messengers' spatial distribution. Continuous thermal agitation induces random events such as macromolecule allosteric transitions, messenger diffusion, and molecular collisions. Therefore, biochemical interactions can be described as a set of strongly coupled causal Markov chains, with random variables specifying internal states such as conformational states and concentrations.

The main result of this paper is a formal demonstration that concentrations of messengers in biochemical systems can represent subjective probability distributions and that conformational changes of macromolecules can perform the fundamental operations required by Bayesian inference.

To reach this goal we demonstrate, on one hand, that Bayesian inference on subjective models is equivalent to the computation of some rational function with nonnegative coefficient (RFNC) of inputs related to sensory data, and, on the other hand, that biochemical cascades may also be seen as computing RFNCs of inputs related to the concentration of primary messengers.

A basic example

The simplest subjective model has only one binary variable S ∈ {0, 1} and a single observation K ∈ K. The structure of the model consists of a prior on S, P (S), and a likelihood of the observation K knowing S, P (K|S):

P (S, K) = P (S)P (K|S) (1)
The prior on S is completely specified by the odds:

P ([S = 1]) P ([S = 0]) = a (2)
It is straightforward to see that P ([S = 0]) = 1/(1 + a) and P ([S = 1]) = a/(1 + a).

The likelihood function P (K | S) allows us to define a mapping g from K to R + such that:

g(k) = P ([K = k] | [S = 1]) P ([K = k] | [S = 0]) (3)
For a given observation k, we apply this function to find the input number x = g(k), from which we can compute the corresponding posterior odds:

y = P ([S = 1] | k) P ([S = 0] | k) = ax (4)
Obviously, this single nonnegative number y completely specifies the posterior distribution with

P ([S = 0] | k) = 1/(1 + y) and P ([S = 1] | k) = y/(1 + y).
This extremely basic subjective model illustrates the two main operations required by probabilistic inference: once an observation k is known, it is first transformed into a nonnegative number x that we call the input (we will later show that in the general case the input is a multidimensional vector of nonnegative numbers). The input is then processed by an RFNC, reduced in this basic model to a multiplication by a constant. The output is again a nonnegative number (in the general case, a multidimensional vector of nonnegative numbers), which specifies the posterior distribution.

We will now show how these two successive operations can be implemented by a biochemical system (see Figure 1). We first consider a population of macromolecular receptors sensitive to some external stimulus. Depending on the type of receptor, the stimulus could be, for instance, a flow of photons, a mechanical strength, or the presence of nutrient in the environment. This stimulus is what the cell observes from its current environment, and its actual particular value can be denoted as k. This stimulus induces a conformational change, and thus the proportion of active receptors. The activity of these receptors results in the release of a primary messenger X within the cell. We call x = [X] the concentration of the primary messenger resulting from the stimulus k. The population of receptors plays exactly the role of the mapping g in the subjective model. Although k is not represented per se in the biochemical system, the output of the population of receptors, i.e., the primary messenger concentration x, represents what really matters for the inference, the likelihood ratio g(k).

We next consider a second population of macromolecules M , each with a single receptor site, specific for the primary messenger X. These macromolecules could be in one of two conformations: in the first, the receptor site is free (M = 0) and in the other the receptor binds a messenger (M = 1). Within a small time interval ∆t, the transition probability from M = 1 to M = 0 is a constant α. The transition from M = 0 to M = 1 requires the presence of one messenger molecule in the vicinity of the receptor site; its probability is therefore proportional to the current messenger concentration x:

P 1→0 = P (M t+∆t = 0 | M t = 1, x) = α P 0→1 = P (M t+∆t = 1 | M t = 0, x) = βx (5)
Each macromolecule in the M population switches randomly between state 0 and state 1 according to the transition probabilities defined above. After a while, the probability of finding the macromolecule in a given state converges to the equilibrium, or steady-state, distribution. In this case, the equilibrium distribution (P * ) follows the detailed balance equation:

P * 1 .P 1→0 = P * 0 .P 0→1 (6)
and the resulting equilibrium distribution is:

P * 0 = α/(α + βx) P * 1 = βx/(α + βx) (7)
When in state 1, that is, when the primary messenger is fixed, the macromolecules release a second messenger Y (concentration y = [Y ]) at a constant rate a P . When in state 0, that is, when the receptor site is free, the macromolecules remove the second messenger at a rate proportional to its concentration a R × y, so that at equilibrium we have:

a P × βx α + βx = a R × α α + βx × y (8)
Finally, the resulting second messenger concentration is simply proportional to the primary messenger concentration:

y = a P a R × βx α (9)
This basic example illustrates the kind of equivalence scheme that can be drawn between the subjective and descriptive models. The macromolecular receptors transform a stimulus k into a likelihood ratio represented by the concentration of primary messengers, while the remainder of the biochemical chain computes a posterior ratio represented by the concentration of second messengers. Obviously, this basic model is quite elementary (see a simulation with realistic values in Figure 1). In the next two sections, we will show that the equivalence scheme is much more general.

Bayesian models and RFNCs

In this paper, we consider only subjective models where unknown variables are discrete, i.e., they have a finite, though eventually huge, number of possible values. The unknown variables can be divided into a set of relevant, or searched variables, denoted S (of size n S ) and a set, eventually empty, of intermediate, or free variables, denoted F (of size n F ). We do not make any restriction on the nature, discrete or not, of the set of observations K. Based on these assumptions, the structure of the general subjective model we consider is:

P (S, F, K) = P (S, F ) × P (K | S, F ) (10)
The right side of this equation is the product of a prior distribution P (S, F ) and a likelihood P (K | S, F ).

For the organism, the problem is to compute the posterior distribution over the relevant variable S knowing a given set of observations k (by convention, we will use small caps for variables with known values):

P (S | k) = F P (S, F ) × P (k | S, F ) S F P (S, F ) × P (k | S, F ) (11)
The denominator on the right side of the equation is a normalization constant, which could be hard to compute if the state spaces (S and F ) are very large. Instead, the posterior distribution is completely determined by a finite set of probability ratios:

P ([S = s] | k) P ([S = 0] | k) = F P ([S = s], F ) × P (k | [S = s], F ) F P ([S = 0], F ) × P (k | [S = 0], F ) (12)
where the particular state [S = 0] is used as a reference, or default state. One can also choose a reference value for the free variables. Providing that the prior P ([S = 0], [F = 0]) and the likelihood

P (k | [S = 0], [F = 0]
) for the chosen reference values are not null, the posterior ratios reduce to:

P ([S = s] | k) P ([S = 0] | k) = f P ([S=s], [F =f ]) P ([S=0] [F =0]) . P (k | [S=s], [F =f ]) P (k | [S=0] [F =0]) f P ([S=0], [F =f ]) P ([S=0], [F =0]) . P (k | [S=0], [F =f ]) P (k | [S=0], [F =0]) (13)

From Bayesian models to RFNCs

Once a given observation k is available, it can be first transformed into an input vector x = g(k) of dimension n SF = n S n F with the following mapping:

∀s ∈ {0, ..., n S -1), ∀f ∈ {0, ..., n F -1) : x s,f = P (k | [S = s] [F = f ]) P (k | [S = 0] [F = 0]) (14)
Obviously, we have x 0,0 = 1 and all other components are nonnegative numbers.

Similarly, the prior P (S, F ) is completely specified by a vector of nonnegative constants a such that:

∀s ∈ {0, ..., n S -1), ∀f ∈ {0, ..., n F -1) : a s,f = P ([S = s] [F = f ]) P ([S = 0] [F = 0]) (15)
Finally, the posterior P (S | k) is completely specified by a vector of dimension n S :

∀s ∈ {0, ..., n S -1) :

y s = P ([S = s] | k) P ([S = 0] | k) (16)
Equation ( 13) can be restated as:

∀s : y s = f a s,f × x s,f f a 0,f × x 0,f (17)
Obviously, we have y 0 = 1 and all other components are rational functions with nonnegative coefficients of the components of the input vector. This demonstrates that the exact inference in the Bayesian model is equivalent to the application of a finite set of RFNCs to an observation-dependent vector of nonnegative numbers.

From RFNC to Bayesian models

Reciprocally, we will now show that from any finite set of RFNCs, one can build at least one subjective probabilistic model such that the posterior probability ratios can be computed with these RFNCs.

Let us first prove this result for a single RFNC h. Clearly, in this case, the search variable S is a binary random variable, and the problem is to specify the subjective model such that the posterior odds is equal to h(x) where x is an input vector depending on the set of observations or known variables k.

We define the complexity of an RFNC as the number of operations (addition, multiplication, and division) required to compute the image h(x) from the components of the input vector.

If complexity(h) = 0, h(x) is either a constant or one component of the input vector. There are obvious Bayesian models associated with both cases, as shown in the basic model section.

If complexity(h) ≥ 1, then it is possible to decompose h into a combination of two simpler, (i.e., lower-complexity) RFNCs h 1 and h 2 with either h

(x) = h 1 (x) + h 2 (x), h(x) = h 1 (x).h 2 (x), or h(x) = h 1 (x)/h 2 (x).
Suppose that we can associate h 1 with a probabilistic model P (S 1 , F 1 , K), where the searched variable S 1 is binary, such that:

P ([S 1 = 1] | k) P ([S 1 = 0] | k) = h 1 (x) (18)
Similarly, we suppose that we can associate h 2 with another probabilistic model P (S 2 , F 2 , K) where S 2 is binary such that the posterior ratio is equal to h 2 (x).

We may assemble these two models using a generative metamodel defined by:

P (S, S 1 , S 2 , F 1 , F 2 , K, K 1 , K 2 , Λ, Γ) = P (S)P (S 1 , F 1 , K 1 )P (S 2 , F 2 , K 2 )P (Λ | S, S 1 , S 2 )P (Γ | K, K 1 , K 2 ) ( 19
)
This metamodel uses two binary variables Λ and Γ, which are called coherence variables (see Chapter 8 in [START_REF] Bessière | Bayesian Programming[END_REF]) and define how the two submodels are combined.

The distribution

P (Γ | K, K 1 , K 2 )
is defined as a Dirac distribution, which means that the same set of observations K is shared by both submodels:

P ([Γ = 1] | K, K 1 , K 2 ) = 1 K=K1=K2 (20)
The distribution P (Λ | S, S 1 , S 2 ) is defined as follows:

P ([Λ = 1] | [S = i], [S 1 = j], [S 2 = k]) = q i,j,k (21)
where the coefficients q i,j,k define how S is related to S 1 and S 2 .

We may now consider the model specified by:

P (S, S 1 , S 2 , F 1 , F 2 , K, K 1 , K 2 | [Λ = 1] , [Γ = 1]) (22)
The posterior ratio on S knowing k for this model is equal to:

P ([S = 1] | k, [Λ = 1], [Γ = 1]) P ([S = 0] | k, [Λ = 1], [Γ = 1]) = q 1,0,0 + q 1,1,0 h 1 (x) + q 1,0,1 h 2 (x) + q 1,1,1 h 1 (x)h 2 (x) q 0,0,0 + q 0,1,0 h 1 (x) + q 0,0,1 h 2 (x) + q 0,1,1 h 1 (x)h 2 (x) (23)
and we can choose the coefficients q i,j,k such that the posterior ratio on S is equal to:

1. either the sum h 1 (x) + h 2 (x) with all coefficients null except q 1,1,0 = q 1,0,1 = q 0,0,0 = 1;

2. or the product h 1 (x).h 2 (x) with all coefficients null except q 1,1,1 = q 0,0,0 = 1;

3. or the quotient h 1 (x)/h 2 (x) with all coefficients null except q 1,1,0 = q 0,0,1 = 1.

Finally, using this generative procedure recursively to derive submodels of null complexity, we can construct for any h a Bayesian model such that:

P ([S = 1] |k) P ([S = 0] |k) = h (x) (24)
Consider now a set of n RFNCs. For each RFNC h i (x), we have shown that there exists at least one (generally several) probabilistic model based on a single searched binary variable S i such that the posterior ratio is equal to h i (x). We now include a global discrete variable S that can take n s = n values, which is related to the binary variables S i by: S = i ⇔ S i = 1 and S j =i = 0. Again, among other possible encoding schemes (possibly more compact), this could be specified by a coherence binary variable Ψ such that P ([Ψ = 1] | S, S 1 , ..., S n ) = 1 if and only if when S = i, S j =i = 0 and S i = 1.

Biochemical cascades and RFNCs

According to the Monod-Wyman-Changeux (MWC) model [START_REF] Monod | On the nature of allosteric transitions: A plausible model[END_REF][START_REF] Changeux | Allosteric mechanisms of signal transduction[END_REF], the activity of a macromolecule depends on its tridimensional tertiary or quaternary structure, which can be in a discrete number of states, typically two, named "tensed" and "relaxed" in the original formulation of the model. The transition probability between these allosteric conformations depends on the status (free or not) of the receptor sites, and the affinity for specific messengers depends on the allosteric conformation. Therefore the state of a macromolecule is defined by the couple (Q, R) where 2), and R = (R 1 , ..., R n R ) is a set of n R binary variables specifying the state of the receptor sites (n R = 2 in the example of Figure 2). Each receptor site can bind a specific messenger, so there are also n R potentially different messenger concentrations represented by the vector x = (x 1 , ..., x n R ). The macromolecule can take 2 n M different states (n M = n Q + n R ) and for each of these states, at equilibrium, the probability of leaving the state is equal to the probability of reaching it:

Q = (Q 1 , ..., Q n Q ) is a set of n Q binary variables specifying the allosteric conformation (n Q = 1 in the example of Figure
∀i : 1 ≤ i ≤ 2 n M , P * i × j =i [P i→j ] = j =i P * j × P j→i ( 25
)
where P i is the probability of being in state i at equilibrium and P i→j is the probability of switching from state i to state j. Some of the P i→j are constants, while some others assume the presence of a messenger and are consequently proportional to the concentration x i of this messenger. Alternatively, we can specify the descriptive model as a Markov chain with a transition matrix T :

P t = T × P t-1 =     P 1→1 P 2 n M →1 P i→j P 1→2 n M P 2 n M →2 n M     P t-1 (26)
where equilibrium is defined by the eventually existing fixed points:

P * = T × P * (27)
It should be noted that (i) the coefficient of T are positive or null and (ii) the coefficients in a column sum to one because they represent the transition probabilities from one state to any other state.

From biochemical cascades to RFNC

The goal of this section is to demonstrate that biochemical cascades without feedback may be seen as computing RFNCs. Specifically, we will prove that, when at equilibrium, the concentrations of second messengers are RFNCs of the concentration of primary messengers.

To reach that goal we will first show that the probability that a macromolecule in its active state is an RFNC of the concentrations (x) of its primary messengers, and we will then demonstrate that when a second messenger is both produced and removed by two populations of macromolecules (see Figure 3), its concentration (y) is an RFNC of the concentrations of the primary messengers.

Probability for a macromolecule to be in its active state as an RFNC of x: In Section 6.1 we presented a general demonstration that the stationary distributions of a Markov model are RFNCs of the transition matrix coefficients. As this demonstration is quite technical, the details have been postponed to the Materials and Methods section. In fact, the coefficients of transition matrix T of the descriptive model are either constants or proportional to the primary messenger concentrations. Therefore, the stationary distribution as well as the probability of being in the active states are RFNCs of these concentrations.

Here, we only present the demonstration for the special case of Figure 2 to give a taste of the complete proof.

The system is made of (i) two primary messengers X 1 and X 2 with concentrations x 1 and x 2 , and (ii) a population of macromolecules with two receptor sites R 1 and R 2 and one active state Q. The allosteric conformation of the macromolecule is defined by the triplet of three binary values (R 1 , R 2 , Q).

For simplicity, we suppose that transitions in the conformation space are restrained to switches of a single binary variable at a time, so that the transition matrix T contains 3 × 8 = 24 nonnull elements (this hypothesis is not necessary for the general demonstration of Section 6.1). Each nonnull element of the transition matrix specifies the probability of changing one particular binary variable from any given initial conformational state. For the receptor state variables, the transition from 0 to 1 involves the presence of one specific messenger in the vicinity of the receptor site. For the eight transitions of concern here, we have transition probabilities proportional to the corresponding messenger concentration x i ; for instance, T 000→100 = α 000→100 × x 1 . The remaining 16 nonnull transition probabilities are constants.

As a second simplifying hypothesis, we suppose that the net chemical and energetic balance of any complete cycle is null (this hypothesis is not necessary for the general demonstration of Section 6.1). Consequently, the Markov chain is reversible and the conformational state distribution converges towards an equilibrium distribution P * (R 1 , R 2 , Q), which satisfies the 24 detailed balance equations of which only eight are independent:

∀ (R 1 , R 2 , Q) ∈ (0, 1) 3      P * (R 1 , R 2 , Q)T R1R2Q→ R1R2Q = P * ( R1 ,R 2 , Q)T R1R2Q→R1R2Q P * (R 1 , R 2 , Q)T R1R2Q→R1 R2Q = P * (R 1 , R2 , Q)T R1 R2Q→R1R2Q P * (R 1 , R 2 , Q)T R1R2Q→R1R2 Q = P * (R 1 , R 2 , Q)T R1R2 Q→R1R2Q (28)
To compute the probability of a given state at equilibrium we can follow any path starting from the reference state (0, 0, 0). For instance, following the path in Figure 2: (0, 0, 0) (no messenger bound), (1, 0, 0) (green triangle bound), (1, 1, 0) (both messengers bound), and (1, 1, 1) (both messengers bound and macromolecule active), we find:

P * (1, 1, 1|x 1 , x 2 ) = T 110→111 T 111→110 T 100→110 T 110→100 T 000→100 T 100→000 P * (0, 0, 0|x 1 , x 2 ) ( 29
)
As each transition coefficient is either a constant or a constant multiplied by one of the two concentrations x 1 or x 2 , the probability of each state at equilibrium is an RFNC of the two concentrations. By marginalizing over R 1 and R 2 , we obtain the probability distribution at equilibrium over the activity states:

P * ([Q = 1]|x 1 , x 2 ) = R1R2 [P * (R 1 , R 2 , [Q = 1]|x 1 , x 2 )] (30)
which is also an RFNC.

In real biochemical systems, the proportion of active sites fluctuates around the mean theoretical value because the population size is finite. An example simulation is shown in Figure 2.

Concentration y of a second messenger as an RFNC of x: The second messenger production is the net result of a reaction (or a chain of reactions) catalyzed by the first population of macromolecules having m P catalytic sites, with a chemical equation of the form:

A → Y + B (31)
The precursor A is a substrate (or set of substrates) present in high concentrations either in the environment or produced in large amount by the basic cellular activity. Generally, A is not involved in cell signaling per se. The concentration of the optional metabolite B is supposed to have no incidence on the catalysis kinetics2 . The production rate Φ P is the product of the number m P × P (Q|x) of active macromolecules by the messenger production rate per active catalytic site3 a P :

Φ P (x) = a P × m P × P * (Q|x) (32)
which is an RFNC of x.

For y, the concentration of the second messenger Y , to reach an equilibrium, the cascade must include a removal reaction (see Figure 3). The messenger removal is the net result of a reaction (or a chain of reactions), catalyzed by a second population of macromolecules, with a chemical equation of the form:

C + Y → D (33)
For C and D, the same assumptions can be made as for A and B. However, the fixation of the messenger on the catalytic site of the removal population requires the presence of messenger in the compartment, so that the removal rate Φ R is proportional to the second messenger concentration:

Φ R (x) = a R × m R × P * (Q |x) × y ( 34
)
where m R is the number of catalytic sites of the second population, a R is the removal rate, Q is the active state of the second population and x is the concentration vector of the primary messengers.

P (Q |x) is an RFNC of x.
At equilibrium, the production and removal rate are equal and we have:

y = a P × m P × P * (Q|x) a R × m R × P * (Q |x) (35)
y is an RFNC of the concentrations appearing in x.

From RFNC to biochemical cascades

Reciprocally, in this section we will finally demonstrate that any RFNC can be computed by a theoretical biochemical cascade. By "theoretical" we mean that the kinetic parameters and overall organization of the cascade necessary to perform the computation of a given RFNC have absolutely no warranty of biological existence or even plausibility.

As we already saw in Section 4.2.2, any RFNC h(x) can be decomposed into either the sum, product or quotient of two simpler RFNCs h 1 (x) and h 2 (x). Suppose that for each function, there exists a (theoretical) biochemical cascade resulting in the release of two second messengers Y 1 and Y 2 , and their concentrations at equilibrium are y 1 = h 1 (x) and y 2 = h 2 (x). Then, we can define two new macromolecules, both with two receptor sites specific to Y 1 and Y 2 ; the first macromolecule releasing and the second removing a new messenger Z. The mean catalytic activity of the first macromolecule is a rational function of y 1 and y 2 with the following general parametric form:

Φ P (y 1 , y 2 ) = a P × n × a 0 + a 1 y 1 + a 2 y 2 + a 3 y 1 y 2 (a 0 + b 0 ) + (a 1 + b 1 )y 1 + (a 2 + b 2 )y 2 + (a 3 + b 3 )y 1 y 2 (36)
Similarly, the mean catalytic activity of the removing macromolecule is a rational function of y 1 and y 2 with the following general parametric form:

Φ R (y 1 , y 2 ) = a R × m × a 0 + a 1 y 1 + a 2 y 2 + a 3 y 1 y 2 (a 0 + b 0 ) + (a 1 + b 1 )y 1 + (a 2 + b 2 )y 2 + (a 3 + b 3 )y 1 y 2 (37)
As we saw above, at equilibrium the concentration of the new messenger Z is given by:

z = Φ P (y 1 , y 2 )/Φ R (y 1 , y 2 ) (38)
We will now show that it is possible to find a set of parameters such that either z = y 1 + y 2 , z = y 1 .y 2 , or z = y 1 /y 2 . Some, but not all, parameters can be set to zero. These parameters determine the probability of the macromolecules being in one of the eight possible states: 5. the active state with no bound messenger is determined by parameters a 0 and a 0 ;

6. the active state with messenger Y 1 bound to the receptor site is determined by parameters a 1 and a 1 ;

7. the active state with messenger Y 2 bound to the receptor site is determined by parameters a 2 and a 2 ;

8. the active state with messengers (Y 1 and Y 2 ) bound to their receptor sites is determined by parameters a 3 and a 3 .

We assume that the resting state has a nonnull probability for both macromolecules. Without loss of generality, we can choose b 0 = b 0 = 1. Obviously, if a parameter is nonnull, meaning that the macromolecule can effectively be in the corresponding state, this implies that there exists a path joining this state to the default state along with all intermediary states that have nonnull parameters. For instance, if a 3 > 0, we must have either (b 3 > 0 and b 1 > 0) or (b 3 > 0 and b 2 > 0) or (a 1 > 0 and b 1 > 0) or (a 2 > 0 and b 2 > 0). Fortunately, there are various sets of parameters that fulfill these requirements and that allow us to compute the sum, the product or the quotient. Here is a possible solution:

• for the sum, we can choose

b 1 = b 2 = a 1 = a 2 = a 0 = 1, b 3 = a 0 = a 3 = b 3 = a 1 = a 2 = a 3 = 0, b 1 = b 2 = 4 and a P n/a R m = 1/2, so that z = Φ P (y 1 , y 2 )/Φ R (y 1 , y 2 ) = y 1 + y 2 ;
• For the product (see Figure 3), we can choose b

1 = b 3 = a 3 = a 0 = 1, b 2 = a 0 = a 1 = a 2 = a 1 = a 2 = a 3 = b 2 = 0, b 1 = 2, b 3 = 4 and a P n/a R m = 1/2
, so that z = y 1 .y 2 ;

• For the quotient, we can choose b

1 = b 2 = a 1 = 1, b 3 = a 0 = a 2 = a 3 = b 3 = a 0 = a 1 = a 3 = 0, b 1 = 2, b 2 = a 2 = a P n/a R m = 1/2, so that z = y 1 /y 2 .
Hence, by recursion, any RFNC (and any set of RFNCs) can be implemented by a cascade of theoretical biochemical processes receiving as inputs a set of primary messengers having concentrations equal (or proportional) to the coefficients of the input vector x.

Discussion

We proposed to look at an aspect of cell signaling and biochemistry from a new viewpoint, assuming that these processes may be seen as performing probabilistic inference. This proposition is founded on mathematical equivalences between descriptive probabilistic models of the interactions of populations of macromolecules and messengers on the one hand, and subjective probabilistic models of the interaction with the environment on the other hand.

The main necessary hypotheses may be summarized as follows.

1. We consider a small volume (order of magnitude: 1µm 3 ) of cell cytoplasm in which thermal diffusion ensures homogeneous concentrations of messengers at the considered time scale (in the range 1 ms to 1 s).

2. We consider one or several populations of allosteric macromolecules with a fixed concentration but variable conformations.

3. We consider one or several messengers with a single conformation but varying concentrations.

4. We consider the conformational changes of the macromolecules, which control the concentrations of the messengers and reciprocally, the concentrations of the messengers, which control the conformational changes of the macromolecules.

Assuming these four hypotheses, we have demonstrated that descriptive models of biochemical systems might be seen as performing probabilistic inference of some well-known and interesting subjective models.

Although the validity of these four hypotheses remains to be discussed, the propositions made in this paper open new perspectives for future directions of study, for instance: can we propose subjective interpretations of descriptive models of other types of biochemical interactions? At the same time and space scale? At different time and space scales?

Validity of the hypotheses?

A critical issue for designing descriptive model is to specify the space and time range in which processes are analyzed.

In the spatial domain, molecular assemblies of nanometer size govern the behavior of organisms that are several meters long. In this paper, we have restrained our analysis to small compartments of a single cell, typically a portion of a dendrite or axon. The order of magnitude of the size of such a compartment is 1µm 3 = 10 -18 m 3 .

Biological events range from about 10 -12 s [START_REF] Knapp | Allosteric action in real time: Time-resolved crystallographic studies of a cooperative dimeric hemoglobin[END_REF][START_REF] Elbert | A milestoning study of the kinetics of an allosteric transition: Atomically detailed simulations of deoxy scapharca hemoglobin[END_REF] for the fastest observed conformational changes to millions of years for gene evolution. We have restrained our analysis to time windows of 1 ms to 1 s, which is long enough to insure homogeneous concentrations of messengers in the considered compartments.

Our description of biochemical system is still too schematic for several reasons.

• The distinction between macromolecules and diffusible messengers is not sharp. Macromolecules also diffuse within membranes [START_REF] Triller | Surface tracking of receptors between synaptic and extrasynaptic membranes: and yet they surface tracking of receptors between synaptic and extrasynaptic membranes: and yet they do move[END_REF] and their spatial distribution, for instance, between cytoplasmic and nuclear compartments, is also controlled by specific biochemical mechanisms [START_REF] Stipanovich | A phosphatase cascade by which rewarding stimuli control nucleosomal response[END_REF]. However, at the space and time scale we are considering, this mobility of macromolecules may be neglected as a first approximation.

• A number of macromolecules directly act on other macromolecules without intermediate messengers.

A clear example in cell signaling is given by the interplay of kinases and phosphatases. This was also discarded from our study for the moment.

• Detailed models of biochemical processes should also take into account the particular geometry of the cell. This includes (3D) phenomena arising in homogeneous volumes (such as the concentration of diffusible messengers), (2D) phenomena arising on a membrane area (such as the distribution of channels), and (1D) phenomena arising mainly along a symmetry axis (such as the propagation of potential along dendrite or axon branches). At present, we have only considered homogeneous concentrations of messengers and fixed distributions of channels, and we have neglected the effect of membrane potentials.

• We have only considered biochemical systems at equilibrium. Clearly, the temporal evolution of concentrations and macromolecular configurations must be further developed, in particular to account for slow reaction and diffusion processes, and for the various roles of biochemical feedback pathways. These dynamic processes could be related to time-evolving probabilistic reasoning, for instance hidden Markov models and Bayesian filters. However, the wide diversity of time scales encountered in biological systems contrasts with the somewhat schematic and oversimplified view of time representation in the usual subjective models. In future work, promising ideas could be the development of a more complex temporal hierarchy, and a more subtle view of the respective roles of memory and temporal reasoning in subjective models.

Subjective interpretation of other biochemical interactions?

Obviously, the above detailed hypotheses are too restrictive. Furthermore, in this work, we have only considered a few possible biochemical interactions among the huge variety of possible interactions. An important task in the near future will be to relax these hypotheses and look for subjective models that could be associated with other kinds of biochemical interactions. Some instances of such perspectives are discussed in the sequel to this section.

At the same time and space scale

Macromolecules with more allosteric states: The allosteric theory [START_REF] Monod | On the nature of allosteric transitions: A plausible model[END_REF], initially developed to account for regulatory enzyme kinetics, postulated that proteins undergo fast, reversible transitions between a discrete number of conformational states. Transitions occur spontaneously, but some are favored by fixation of ligand to specific receptor sites. This model has been successfully applied to a large number of fundamental macromolecular assemblies, such as hemoglobin, ionic channels, and nuclear or membrane receptors (see [START_REF] Changeux | Allosteric mechanisms of signal transduction[END_REF] for a review).

The number of conformational states and the variety of controlled mechanisms for conformational changes can be relatively high. As an example, DARPP-32, a key macromolecule in the integration of dopamine and glutamate inputs to the striatal GABAergic neurons, exhibits four phosphorylation sites, thus 16 conformational states, and each phosphorylation is controlled by a different chemical messenger pathway [START_REF] Fernandez | DARPP-32 is a robust integrator of dopamine and glutamate signals[END_REF]. Inactivation of rhodopsin in vertebrate photoreceptors involves up to 12 phosphorylation sites, which are all controlled by the intracellular calcium concentration. In [START_REF] Houillon | The probabilistic cell: Implementation of a probabilistic inference by the biochemical mechanisms of phototransduction[END_REF] we described a first application of our approach to the vertebrate phototransduction biochemical cascade. Augmenting the number of conformational states opens very exciting perspectives on the complexity of the computation that a single population of macromolecules could perform, but the corresponding subjective models are still to be proposed.

Macromolecules with several receptor sites for the same messenger: Macromolecules with several receptor sites for the same messenger are very frequent in biochemistry. A number of allosteric macromolecules, including ionic channels, are composed of several subunits, and some of them are identical. The presence of a receptor site for a specific messenger on each subunit makes the whole macromolecule controlled by the second, third, or fourth power of the concentration. As a consequence, the catalytic activity is a highly nonlinear function of the messenger concentration, and may exhibit very sharp sensitivity.

In terms of subjective models, this means that similar observations, converted into likelihood ratios, are performed several times to infer the posterior distribution. This might be a simple and elegant way to enrich the computational complexity performed by macromolecules without requiring high dynamic ranges of messenger concentration. In the present work, we have assumed a simple proportional relationship between likelihood ratios and messenger concentrations. The presence of multiple receptor sites for the same messenger, and more generally the multimeric structure of some macromolecules can be interpreted as a nonlinear coding of the likelihood ratio, which could be better adapted to the biological constraints. For instance, a tenfold increase in messenger concentration could correspond to a likelihood ratio multiplied by 10 4 for a tetrameric receptor.

Allosteric changes governed by other events instead of chemical messengers: Allosteric changes may be caused by other types of events besides chemical messengers such as electrical or mechanical ones. These events are not considered in the present model and should be studied in future.

Single-channel voltage clamp recordings [START_REF] Neher | Single-channel currents recorded from membrane of denervated frog muscle fibres[END_REF][START_REF] Sakmann | Single-Channel Recording[END_REF] have revealed several important characteristics of channels: (i) currents through isolated channels, and thus channel conductance, alternate between discrete values; (ii) transitions between current/conductance values are random brisk events; (iii) the transition probabilities can be modulated by pharmacological and biological agents (neurotransmitters, second messengers), ions (calcium), or membrane potentials. In agreement with the allosteric theory, the current descriptive model of ion channels [START_REF] Colquhoun | On the stochastic properties of single ion channels[END_REF][START_REF] Laüger | Conformational transitions of ionic channels[END_REF][START_REF] Chung | Coupled markov chain model: characterization of membrane channel currents with multiple conductance sublevels as partially coupled elementary pores[END_REF] is that of a finite state Markov model, similar to the one we used in this paper. Some transitions depend on the presence of a specific messenger in the vicinity of receptor sites, either in the extracellular domain (ionotropic receptor-channels) or in the intracellular domain, where receptor sites are specific to second messengers or calcium transporters. In voltage-dependent channels, transition probabilities are controlled by the membrane potential.

Local cascades and feedbacks: Metabotropic receptors are transmembrane macromolecular assemblies with a receptor site in the extracellular domain and an activatable site in the intracellular domain. Chemical binding of the neurotransmitter with the extracellular receptor site induces allosteric conformational change of the macromolecule, which activates the intracellular site (an activatable enzyme site such as protein kinase, protein phosphatase, or G-protein release). More generally, macromolecules can be activated by various events like mechanical strength, photoisomerization, pheromones, or odor detection. The intracellular activity results in the release of a primary messenger (e.g., the G α subunit in the G-protein-dependent signaling), which triggers the production or elimination of a diffusible second messenger (cyclic nucleotides, inositol triphosphate, . . . ), which in turn acts on ionic channels.

In the complex molecular chain from the primary receptor to the set of ionic channels, several allosteric macromolecules intervene. Most of them have receptor sites for calcium, or are calcium transporters like calmodulin, thus allowing feedback regulation. As discussed above, the role of feedback should be understood in terms of the dynamics of the systems, including long-term adaptation.

At different time and space scales

Different dynamics to store information: In the core of our model, we postulate that the existence of well-separated fast rates (10 4 to 10 5 s -1 ) and relatively slow rates (10 to 10 3 s -1 ) of biochemical reactions is fundamental. It allows the Markov processes of molecular collision and configuration changes generated by thermal agitation to converge to quasistationary states that approximate probabilistic inference of subjective models. We have not yet considered the very slow dynamics of some allosteric changes (like desensitization) or the diffusion of large molecules in membranes, among many other slow biochemical events. These mechanisms can clearly be used for accumulating and storing information over long periods, which is a key computational capacity for adaptation and learning.

The role of membrane potential at cell scale: The overall effect of all ionic fluxes can be summed up in the membrane potential kinetics. Though it is rather unusual to include membrane potentials in the set of messengers, it seems appropriate because the kinetic equations of membrane potential are similar to those for chemical messenger concentrations. Moreover, the membrane potential controls macromolecule conformational transition in a similar way, although constrained to a particular, but important, class of macromolecules, namely the voltage-gated channels [START_REF] Bezanilla | The voltage sensor in voltage-dependent ion channels[END_REF]. In turn, ionic channels change the membrane potential similarly to the way that activatable enzymes control chemical messenger concentration. Membrane potentials propagate much faster than messengers diffuse and can transmit the result of a given local computation to the whole cell. Consequently, including membrane potentials into our proposed framework would extend the space scale from a compartment to the whole cell.

Unicellular organisms: At this molecular description level, our proposal applies not only to braincontrolled complex organisms and not only to small neuron networks, but also to unicellular organisms. Simple organisms like Paramecium or Euglena gracilis have limited numbers of sensors, and a greatly reduced repertoire of actions, but they must nevertheless adapt their behavior to an even more unpredictable environment. The efficiency of probabilistic reasoning with an incomplete model of the world applies equally to these simple organisms. Furthermore, the biochemical mechanisms that we propose for implementing probabilistic computation are already effective in controlling the behavior of eukaryotes. Among many other examples, it has been shown [START_REF] Iseki | A blue-light-activated adenylyl cyclase mediates photoavoidance in euglena gracilis[END_REF] that the photoavoidance behavior of the microalgae (E. gracilis) is mediated by concentration changes of cyclic adenosine monophosphate, a second messenger known to be involved in olfaction and many neuronal signaling pathways of multicellular animals.

Information transmission in multicellular organisms: Excitable cells, particularly neurons, differentiate from other cells in complex multicellular organisms, from cnidarians to all bilaterians. These excitable cells provide a distant and discrete mode of messenger flux control: both chemical or ionic diffusion and passive electrical propagation become very ineffective at long distances. Note, however, that slow, but distant, signal propagation through active biochemical processes without action potentials has been recently discovered [START_REF] Fasano | Neuronal conduction of excitation without action potentials based on ceramide production[END_REF]. On the contrary, active regeneration of action potential thanks to voltagegated channels allows the message to be transmitted unchanged at high speed along axonal branches up to the presynaptic terminals, where it is converted into chemical signals. Distant interactions between macromolecular assemblies by action potential propagation allow multicellular organisms to reach sizes and speeds far greater than the limit imposed by passive diffusion process. This constitutes an obvious gain for long-distance communication. However as opposed to membrane potential and messengers' concentrations which are graded and local signals, the action potentials are all-or-none signals. Therefore long distance communication is achieved at the cost of signal impoverishment. According to our view, fast-spike propagation must be completed by local and graded signal processing involving complex biochemical interactions that constitute a fundamental mechanism for probability computation.

6 Materials and Methods We consider a set {X i } 1≤i≤n of states with the transition probabilities among those states T . The dynamics for the Markov chain with initial occupancy P 0 is thus P t = T P t-1 :

P t =        T 11 • • • T 1n . . . . . . T ij . . . T n1 • • • T nn        P t-1
where the probability of a transition from X j to X i is T ij . If the Markov process is finite and irreducible, then according to the Perron-Frobenius theorem, it has a unique stationary distribution π [START_REF] Stroock | An Introduction to Markov Processes[END_REF]:

π i = det((I -T ) {i} ) n k=1 det((I -T ) {k} ) , (39)
where M {k} = (I -T ) {k} denotes the matrix obtained after removing the term m ii and its corresponding column and row. For every X j , we have n i=1

T ij = 1. We note that s j = 1 -T jj is positive as

s j = n i=1 i =j T ij . M {k} =    s i1 • • • -T i1jn-1 . . . . . . . . . -T in-1j1 • • • s in-1   
By the Leibniz formula:

detM {k} = σ∈S (σ) n j=1 j =k
M σ(j)j , where S are the permutations over indexes {1, • • • , n}\k. Because of the form of the diagonal coefficients, this can be rewritten as:

detM {k} = a∈A {k} λ(a) n j=1 j =k
T a(j)j , where A k is the set of applications:

a : {1, • • • n}\k → {1, • • • n} j → a(j) = j
, that is, the set of applications on indexes from 1 to n except k, leaving no indexes invariant and with k as a possible image. We consider elementary matrices a M {k} with ∀j ∈ {1,

• • • n}\k ∀i ∈ {1, • • • n}\k s.t. i = a(j) a M {k} ij = 0.
The determinant of such a matrix is reduced to:

det a M {k} = λ(a) n j=1 j =k
T a(j)j . Note that the coefficient λ(a) is the same for det a M {k} and detM {k} , so that: detM {k} = a∈A {k} det a M {k} .

Almost-diagonal matrices

For a matrix (hereafter referred as an almost-diagonal matrix) having positive coefficients on the diagonal and at most one other single nonzero coefficient in each column, this additional coefficient being opposite to the diagonal term, we denote by K(M ) the number of columns having a single nonzero coefficient. The elementary matrices defined above are almost-diagonal with K( a M {k} ) = #{j : a(j) = k} and applications a M {k} can then be studied based on their K values.

If K = 0 s j = T a(j)j and the matrix is reduced to:

a M {k} =                T a(1)1 • • • • • • 0 • • • 0 0 -T a(j)j . . . . . . . . . . . . -T a(n)n -T a(1)1 . . . 0 0 T a(j)j . . . . . . . . . . . . 0 0 • • • 0 • • • T a(n)n               
For such a matrix, the terms in each column sum to zero so that det a M {k} = 0 and then λ(a) = 0.

If K = n -1 In this case, the matrix is diagonal and

det a M {k} = n j=1 j =k T kj so that λ(a) = 1. If K ∈ {1, • • • , n -2}
There is at least one column with only its diagonal term nonzero. After Cramer expansion along one such column (for example the smallest j such that a(j) = k), det a M {k} = T kj det a M {k}{j} . The new matrix a M {k}{j} is of dimension n -2 and is almost diagonal. We name K 1 ( a M {k} ) ∈ 0, • • • , n -2 the number of columns having a single nonzero coefficient (on the diagonal).

If ∃i|a(i) = j, K 1 ( a M {k} ) = K( a M {k} ), otherwise K 1 ( a M {k} ) < K( a M {k} ), we then iterate this operation generating the sequence (K s ( a M {k} )) 1≤s≤n-2 . As the K s are bounded by the matrix dimension ∀s ∈ {1, • • • , n -2}, K s ( a M {k} ) ≤ ns -1, there is either an s = s + for which K s + (a) = ns -1 and then λ(a) = 1 (as shown in the case K = n -1) or an s = s -for which K s -(a) = 0 and then λ(a) = 0 (as shown in the case K = 0).

6.1.3

Stationary distributions of Markov models are RFNCs of the transition matrix coefficients.

We showed that the determinants of M {k} = (I-T ) {k} can be decomposed as detM {k} = a∈A {k} det a M {k} and ∀a ∈ A {k} , λ(a) = 1 or λ(a) = 0, so these determinants are polynomials in the transition probabilities with nonnegative coefficients. It then follows from (39) that the stationary distributions for a Markov chain are RFNCs of the transition probabilities.

Simulations

Algorithms

We performed simulations with a software package written in C++. To simulate the descriptive models, we used the direct version of the exact stochastic Gillespie algorithm [START_REF] Gillespie | Exact stochastic simulation of coupled chemical reactions[END_REF][START_REF] Gillespie | Stochastic simulation of chemical kinetics[END_REF]. The current state of the biochemical system is a vector of natural integers X t . Each component X t i corresponds to the current number of molecules in a particular conformation in a particular compartment. The reactions are considered as random events R t . Each reaction/event can occur with an instantaneous rate (or propensity) proportional to the number of reactants, e.g., if i and j are the index of the reactants involved in reaction r:

a(r, t) = c(r)X t i X t j (40)
The next reaction is drawn randomly from the histogram of reaction rates:

r t ∼ P (R t = r) = a(r, t) a 0 (t) , with a 0 (t) = r a(r, t) (41)
The elapsed time at which the next reaction occurs after time t is drawn from the exponential distribution with intensity equal to the sum of reaction rates: ∆t ∼ a 0 (t) -1 e -a0(t)∆t (42)

The state vector is then updated by the stoichiometric vector of the reaction:

X t+∆t = X t + S(r t ) (43)
The components of the stoichiometric vector are integers defining the change in reactant number resulting from the occurrence of a given reaction. For instance, if the reaction r involves the collision of reactants i and j, which results in the formation of a new reactant k, then we have S i (r) = S j (r) = -1 and S k (r) = +1. In addition, some components of the state vector are set to prespecified values. They constitute the inputs to the biochemical process.

Bayesian models including Bayesian filters are simulated by computing the exact inference, i.e., by applying the Bayesian rule and marginalization rules on histogram distributions. The integer numbers of molecules of each species are computed for a diffusion volume of 1 µm 3 . The acetylcholine concentration is fixed to 2 µM (i.e., n 1200 molecules) for the first 10 ms then jumps to 20 µM for the next 5 ms, then drops back to the initial value.

Parameters

For multiple realizations of the simulation, a differential equation for the mean fraction of activated macromolecules can be derived from the master equation (see [START_REF] Gardner | Handbook of stochastic methods[END_REF]):

dq(t) dt = -αq(t) + βx(1 -q(t)). (44)
Similarly, the dynamics of the variance is: ). There are nine molecular species (the ACh messenger plus the eight conformational states of the channel-receptor), and 24 reactions. Both receptor sites have the same fixation (β) and removal rates (α 0 and α 1 ). The transition rates between Q 0 and Q 1 states depend on the number i (0, 1 or 2) of fixed ligands. We simulate the transitions of 200 macromolecules using the Gillespie algorithm. Half of the molecules are involved in the production of Y and may be activated when both receptor sites are occupied and receptor binding is sequential (first X 1 , then X 2 ). The other half are involved in the degradation of Y and may be activated only when none of the receptor sites is occupied. Messenger fixation also occurs in a sequential manner. Fixation rates are all assumed to have the same value β and removal rates have the same value α P,R 1,2 , except for the removal of X 2 in the producing macromolecule α R 2 . Activation and deactivation rates also have the same values k a P,R and k d P,R except for the deactivation rate of the removing macromolecule k d R . The dynamics for the production and removal of Y are simulated through integration of the following differential equation:

dv(t) dt = α(q(t) -2v(t)) + βx(1 -q(t) -2v(t)) (45)
τ dy(t) dt = -Φ R (x 1 , x 2 )y(t) + Φ P (x 1 , x 2 ). ( 46
)
7 Acknowledgments

This work is part of the European FP7 project BAMBI (Bottom-up Approaches to Machines dedicated to Bayesian Inferences). The authors warmly thank Emmanuel Mazer, Julie Grollier, and Damien Qurelioz for their fruitful comments. The gray line shows the fraction of bound molecules as simulated with the Gillespie algorithm. The green dots show the ensemble average (± one standard deviation) of these trajectories over 500 repetitions. The orange curve shows the same variable as obtained by the continuous time dynamics of the mean and variance (see Section 6.2.2 for the derivation of these equations). Parameters are chosen according to the acetylcholine receptor described in [START_REF] Edelstein | A kinetic mechanism for nicotinic acetylcholine receptors based on multiple allosteric transitions[END_REF] (see Table 1).

1 . 1 ; 3 .

 113 the resting (default) state with no bound messenger and no catalytic activity is determined by parameters b 0 and b 0 ; 2. the inactive state with messenger Y 1 bound to the receptor site is determined by parameters b 1 and b the inactive state with messenger Y 2 bound to the receptor site is determined by parameters b 2 and b 2 ; 4. the inactive state with both messengers (Y 1 and Y 2 ) bound to their receptor sites is determined by parameters b 3 and b 3 ;

6. 1 6 . 1 . 1

 1611 Stationary distributions of Markov models are RFNCs of the transition matrix coefficients: a formal demonstration Decomposition for the stationary distribution of a Markov model.

Figure 1 :

 1 Figure 1: We simulate a population of N = 100 receptor sites (R) specific for the messenger acetylcholine (X) and the following two reactions: Fixation: R + X → RX Rate coefficient: β Removal: RX → R + X Rate coefficient: αThe integer numbers of molecules of each species are computed for a diffusion volume of 1 µm 3 . The acetylcholine concentration is fixed to 2 µM (i.e., n 1200 molecules) for the first 10 ms then jumps to 20 µM for the next 5 ms, then drops back to the initial value.For multiple realizations of the simulation, a differential equation for the mean fraction of activated macromolecules can be derived from the master equation (see[START_REF] Gardner | Handbook of stochastic methods[END_REF]):

Figure 2 :

 2 Figure 2: We simulate a population of m = 100 NACh channel-receptors with two receptor sites and one catalytic siteQ 0 and Q 1 ([Q = 0] and [Q = 1]). There are nine molecular species (the ACh messenger plus the eight conformational states of the channel-receptor), and 24 reactions. Both receptor sites have the same fixation (β) and removal rates (α 0 and α 1 ). The transition rates between Q 0 and Q 1 states depend on the number i (0, 1 or 2) of fixed ligands.

Figure 3 :

 3 Figure3: We simulate the transitions of 200 macromolecules using the Gillespie algorithm. Half of the molecules are involved in the production of Y and may be activated when both receptor sites are occupied and receptor binding is sequential (first X 1 , then X 2 ). The other half are involved in the degradation of Y and may be activated only when none of the receptor sites is occupied. Messenger fixation also occurs in a sequential manner. Fixation rates are all assumed to have the same value β and removal rates have the same value α P,R 1,2 , except for the removal of X 2 in the producing macromolecule α R 2 . Activation and deactivation rates also have the same values k a P,R and k d P,R except for the deactivation rate of the removing macromolecule k d R . The dynamics for the production and removal of Y are simulated through integration of the following differential equation:

m: size of the population m : 100 NFigure 1 .

 1001 Figure 1. (Left) Schema of the reaction involved in a single receptor site. The two constants α, β are the kinetic constants for the release and fixation of the messenger, respectively. (Right) Dynamics of receptor binding when the ligand concentration is increased for 5 ms for a population of 100 receptors.The gray line shows the fraction of bound molecules as simulated with the Gillespie algorithm. The green dots show the ensemble average (± one standard deviation) of these trajectories over 500 repetitions. The orange curve shows the same variable as obtained by the continuous time dynamics of the mean and variance (see Section 6.2.2 for the derivation of these equations). Parameters are chosen according to the acetylcholine receptor described in[START_REF] Edelstein | A kinetic mechanism for nicotinic acetylcholine receptors based on multiple allosteric transitions[END_REF] (see Table1).

Figure 2 .Figure 3 .

 23 Figure 2. (Left) Schema of a macromolecule with two receptor sites (red circle and green triangle) and one catalytic site (yellow). (Right) Simulation of this macromolecule when the concentrations of both ligands are increased from 2µM to 20µM during 5 ms. The vertical axis shows the fraction of activated catalytic sites. The same conventions as in Figure 1 are used for the curves. The black line shows the activated fraction of catalytic sites in the stationary state. See Table3for the definition of parameters.

Table 1 .

 1 Parameter values for Figure 1 (Ligand binding.)

	Parameter name	Value
	β: fixation rate	β : 150 µM -1 s -1
	α: release rate	α : 8000s -1
	x down , x up : ligand concentrations	x down : 2µM , x up : 20µM
	m: size of the population	m : 100
	N : Number of realizations	N : 500
	dt: time step for ODE integration	dt : 0.005ms

Table 2 .

 2 Parameter values for Figure2(acetylcholine receptor.) 0 , α 1 : removal rates α 0 : 8000.s -1 , α 1 : 8.64 s -1 k a m : activation rate when m ligands are bound k a 0 : 0.54, k a 1 : 130, k a 2 : 30000s -1 k d m : deactivation rate when m ligands are bound k d 0 : 10800, k d 1 : 2808, k d 2 : 700s -1 x down , x up : ligand concentrations x

	Parameter name	Value
	β: fixation rate	β : 150µM -1 .s -1

α down : 2µM , x up : 20µM m: size of the population m : 100 N : Number of realizations N : 500 dt: time step for ODE integration dt : 0.005ms

Table 3 .

 3 Parameter values for Figure3(Abstract reaction for product computation.)

	Parameter name	Value
	β: fixation rate for
	ligands	β = 1000µM -1 .s -1
	α R i : release rates for
	ligand x i (i ∈ 1, 2) on the
	macromolecule involved in production of Y	β P 1 : 1000µs -1 , β P 2 : 2000µs -1
	α R : removal rate for
	ligands	α R : 1000µs -1
	k a P : activation rate for the
	macromolecule involved in production of Y	k a P : 2000s 1
	k a R : activation rate for the
	macromolecule involved in degradation of Y	k a R : 2000s 1
	k

d : deactivation rate k d : 1000s -1 x 1 , x 2 : ligand concentrations x i ∈ [1µM, 10µM ]

At present, Bernstein Center for Computational Neuroscience, Berlin, Germany

However, there are several well-known exceptions to this rule. For instance, the cleavage of phospholipid PIP2 by the enzyme phospholipase C leads to the release of two messengers: diacylglycerol and inositol triphosphate that in turn target different biochemical systems.

We assume that, at the temporal scale considered here, the catalytic activity per active site is constant


Pascal Moyal
email: pascal.moyal@ird.fr

Is Acrapex fletcheri Laporte really an Acrapex (Lepidoptera : Noctuidae)?

This study shows that the species named Acrapex fletcheri Laporte does not belong to the genus Acrapex. This is clear from the adult morphology and the shape of the cremaster of the pupa. This is moreover confirmed by a Bayesian molecular phylogenetic analysis from two mitochondrial genes. In this analysis, A. fletcheri is the sister species of Sesamia sylvata Janse, a species that was proved in a previous study not to belong to the genus Sesamia but to a new genus that has to be described. Both species have similar habitus, but the deepness of the node of their branching as well as the different cremaster suggest they belong to two different genera, and that A. fletcheri should become the type species of another new genus. The genus Acrapex, represented by Acrapex syscia Fletcher in the analysis, is strongly connected to the genera Conicofrontia and Poecopa, which indicates common ancestry that is difficult to deduce from the rather different adult morphology of these species. However, these three genera have in common the absence of cremaster in the pupa. Another grouping, between the genera Buakea and Manga, is in accordance with the shape of the cremasters, but it is little supported (posterior probability=0.68). This phylogeny suggests that the ancestor of the African Sesamiina had likely a cremaster with two spines.

Introduction

The study on the cremasters of the pupae of Sesamiina [START_REF] Moyal | African Sesamiina (Lepidoptera: Noctuidae): the cremaster and the genus[END_REF] showed that, except for Pirateolea, this structure was characteristic of the genus. The pupa of Acrapex fletcheri Laporte has a peculiar bifid cremaster (figure 1), different from that of the other species of Sesamiina, including other Acrapex, such as Acrapex syscia Fletcher and the close species Acrapex brunnea Hampson, who have no cremaster [START_REF] Moyal | African Sesamiina (Lepidoptera: Noctuidae): the cremaster and the genus[END_REF]. This suggests that A. fletcheri might not belong to the genus Acrapex or to any other known genus of Sesamiina.

This study aims at presenting the genus Acrapex and at checking, using morphological and molecular data, if A. fletcheri is really an Acrapex.

Materials and Methods

Larvae of the studied stem borers were collected within their host plants in East Africa and southern Africa. They were reared on artificial diet (2) until pupation and the species were identified after emergence of the adults. The specimens of A. fletcheri were collected in the same place as the type described by Laporte (3), Mau escarpment, in Kenya. The pictures or slides of the types were obtained from museums (London, Paris, Pretoria).

Sequence data were obtained from fragments of two mitochondrial genes, Cytochrome b (Cytb) (915 nucleotides), Cytochrome c Oxidase, subunit 1 (CO1) (892 nucleotides), Total DNA was extracted using a Qiagen DNeasy tissue kit (Qiagen GmbH, Germany).The PCR rofile for the genes was : initial denaturation for 5 min at 92 °C; 35 cycles of denaturation for 1 min at 92 °C, annealing for 1.30 min at 46 °C, extension for 1.30 min at 72 °C; and final extension for 5 min at 72 °C. The reaction mixture contained 3 mM MgCl2, 0.4 µM primers, 0.24 µM dNTPs, 2 U of Promega Taq polymerase and 100 ng of DNA per 50 µl of reaction mixture. The primers used were : for cytb, CP1 (5'-GATGATGAAATTTTGGATC-3') (modified from (4)) and TRs (5'-TCT ATCTTATGTTTTCAAAAG-3') [START_REF] Janse | The Moths of South Africa[END_REF] ; for CO1, Ron (5'-GGATCACCTGATATAGCATTCCC-3') and Hobbes (5'-AAATGTTGNGGRAAAAATGTTA-3') [START_REF] Monteiro | Phylogeny of Bicyclus (Lepidoptera: Nymphalidae) inferred from COI, COII, and EF--1alpha gene sequences[END_REF]. The PCR product was then purified using the Qiagen QIAquick PCR purification kit (Qiagen GmbH, Germany). Sequencing reactions were carried out using the Sanger dideoxy method (6), and, finally, sequences were run and detected on an ABI 377 automated séquencer.

The sequences were aligned with Multalin software [START_REF] Corpet | Multiple sequence alignment with hierarchical clustering[END_REF]. The nucleotide substitution model used was GTR with proportion of invariable sites and gamma distribution. Phylogenies were inferred using a Bayesian method carried out using MrBayes 3.1 [START_REF] Ronquist | MRBAYES 3: Bayesian phylogenetic inference under mixed models[END_REF].

Results and discussion

1-The genus Acrapex

The genus Acrapex was created by Hampson in 1894. The type species is Acrapex prisca, from India (figure 2). In his major work of 1910 (9), Hampson described six new species of Acrapex and the other nine species already included in the genus, and defined again the genus. In particular, he mentioned two main features: « antennae of male typically ciliated » and « Forewing with the apex produced and acute, the termen obliquely curved and not crenulated ». He gave an identification key based on the male antennae (biserrate and fasciculate, with quadrate uniserrations, annulated and ciliated). Only the adult habitus was used at the time to classify the species.

In 1939, Janse [START_REF] Janse | The Moths of South Africa[END_REF] gave a new definition of the genus based on the species Acrapex aenigma Felder & Rogenhofer, of which he had specimens (figure 3). He described precisely the habitus, and, for the first time in the genus, the male genitalia. He told that 17 species had already been described in the genus, among which 14 are African and eight found in South Africa. He described five new species. Therefore, at this time, the species number in the genus amounted to 22. In 1967, Viette [START_REF] Viette | Insectes Lépidoptères Noctuidae Amphipyrinae (Part.) et Melicleptriinae. Faune de Madagascar[END_REF] reported four species in Madagascar, among which two were new. He gave a full description of the genus, which included the habitus and the male and female genitalia. According to him, the antennae are typically ciliated in the males and the forewings have a pointed apex and an exterior margin regularly curved.

The new species described by Viette are very close to the species typical of the genus, which were described by Hampson and Janse. Likewise, his definition of the genus genitalia is that of the classical species. Nevertheless, Janse [START_REF] Janse | The Moths of South Africa[END_REF] had placed previously in the genus several species whose genitalia were noticeably different from Viette's description. He mentioned it in his description of Acrapex mystica, where he told that the genitalia in the genus Acrapex fall into three groups.

Berio [START_REF] Berio | Nuove species e generi di noctuidae africane e asiatische e note sinonimiche. Parte II[END_REF] was faced with this diversity when he studied 16 species he placed in the genus Acrapex. He explained he had no possibility to observe the types and that he was not absolutely certain that some of the species he described were different from species that were previously described. Their genitalia were close to those of species described by Janse, but their habitus was very different. He decided therefore to describe these species as new ones until someone can check whether they are synonymous or not. Moreover, he added a fourth group to the three groups of genitalia described by Janse, because none of the pictures of the latter seemed to correspond to these genitalia. Following Berio, Laporte (3) placed many new species in the genus Acrapex, among which several would be misplaced according to Holloway [START_REF] Ronquist | MRBAYES 3: Bayesian phylogenetic inference under mixed models[END_REF]. In my study on the genera Sciomesa and Carelis [START_REF] Moyal | Systematics and molecular phylogeny of two African stem borer genera, Sciomesa Tams & Bowden and Carelis Bowden (Lepidoptera: Noctuidae)[END_REF], I transferred one of these species in Sciomesa (Acrapex boulardi Laporte became Sciomesa boulardi) and I identified four cases of synonymy: Acrapex congitae Laporte synonymous with Sciomesa boulardi; Acrapex sparsipuncta Laporte synonymous with Sciomesa excelsa (Laporte); Acrapex bryae Laporte and Acrapex fayei Laporte synonymous with Feraxinia jemjemensis (Laporte).

The genus Acrapex appears then still confused: some of the species described by Berio may be synonymous with species described by Janse; the genitalia diversity raises the question of the limits of the genus; many species considered as different by Laporte might be synonymous. The revision of the genus using molecular data would be of great interest. I will limit myself hereinafter to the case of Acrapex fletcheri Laporte.

2-Acrapex fletcheri Laporte

The species type was collected in Kenya, in Mau-escarpment, in 1975 (3). It is shown in figure 4. As it is in rather bad condition, I show, in figure 5, a specimen obtained from a larva collected in its host plant in the same place. This species is more robust and stockier than the typical Acrapex such as those shown in figures 1 and 2. Its forewings are less oblong. Moreover, its antennae are strongly bipectinate, whereas, according to Hampson, Janse and Viette, the antennae of Acrapex are either ciliated without serration or serrate. Its male genitalia are also different from those of the other Acrapex, with the presence of a clasper, and a long and thin cucullus versus a broad, short and rounded one in the other species.

The morphological features of the adult, as well as the cremaster of the pupa, suggest therefore that this species does not belong to the genus Acrapex.

The phylogenetic tree that resulted from the analysis of the molecular data is shown in figure 6. This tree shows that the species called Acrapex fletcheri has no recent common ancestor with the species typical of the genus, Acrapex syscia, but is the sister of a species misplaced in the genus Sesamia, S. sylvata. This group is rather well supported (posterior probability=0.81). The deep node and the long branches as well as the different cremasters suggest that both species belong to two different genera that have to be created. The likeness of the habitus of the adults (figure 7) suggests that both genera should make up a monophyletic genus group. Furthermore, the phylogenetic tree shows that the genus Acrapex, represented by A. syscia, is strongly connected to a group made up of the genera Conicofrontia and Poecopa (posterior probability=0.95). Yet, the morphology of the adults of these species looks rather different, from the frail Acrapex to the robust Conicofrontia. But they share in common the absence of cremaster in the pupa [START_REF] Moyal | African Sesamiina (Lepidoptera: Noctuidae): the cremaster and the genus[END_REF]. Another grouping, between Buakea and Manga, although much less supported (posterior probability=0.68), is in line with the shape of the cremaster.

This phylogeny shows therefore that the species called Acrapex fletcheri does not belong to the genus Acrapex and should become the type species of a new genus. It shows also that most genus groups have at least one genus whose pupa has a cremaster with two spines, which suggests that the ancestor of the African Sesamiina had such cremaster. The absence of cremaster seems to be a more recent and secondary character, which is ancestral in the group Acrapex-Conicofrontia-Poecopa, whereas, regarding the Sciomesa genus group and perhaps Poeonoma serrata, it would result from an evolutionary convergence.

Figure 1 .

 1 Figure 1. Cremaster of the pupa of Acrapex fletcheri Laporte.

Figure 2 .

 2 Figure 2. Acrapex prisca Hampson: type of the species and of the genus Acrapex.

Figure 3 .

 3 Figure 3. Acrapex aenigma Felder & Rogenhofer. A. Male adult. B. Male genitalia.

Figure 5 .

 5 Figure 5. Acrapex fletcheri: male obtained from larva.

Figure 6 .

 6 Figure 6. Phylogenetic tree resulting from the bayesian analysis of the two mitochondrial genes. Numbers indicate Bayesian posterior probabilities. The quotation marks indicate species misplaced in the Sesamia genus.

Figure 7 .

 7 Figure 7. Habitus of the two sister species in the molecular phylogeny. The quotation marks indicate that the species is misplaced in the genus.