H Tembine
email: tembineh@gmail.com

P Vilanova

M Debbah

Noisy mean field game model for malware propagation in opportunistic networks

Keywords: stochastic game, malware propagation, mean field

In this paper we present analytical mean field techniques that can be used to better understand the behavior of malware propagation in opportunistic large networks. We develop a modeling methodology based on stochastic mean field optimal control that is able to capture many aspects of the problem, especially the impact of the control and heterogeneity of the system on the spreading characteristics of malware. The stochastic large process characterizing the evolution of the total number of infected nodes is examined with a noisy mean field limit and compared to a deterministic one. The stochastic nature of the wireless environment make stochastic approaches more realistic for such types of networks. By introducing control strategies, we show that the fraction of infected nodes can be maintained below some threshold. In contrast to most of the existing results on mean field propagation models which focus on deterministic equations, we show that the mean field limit is stochastic if the second moment of the number of object transitions per time slot is unbounded with the size of the system. This allows us to compare one path of the fraction of infected nodes with the stochastic trajectory of its mean field limit. In order to take into account the heterogeneity of opportunistic networks, the analysis is extended to multiple types of nodes. Our numerical results show that the heterogeneity can help to stabilize the system. We verify the results through simulation showing how to obtain useful approximations in the case of very large systems.

Introduction

In modern times, the massive use of information formed an interconnected global society of billions in which communication systems are vital infrastructures. Among them, wireless communications is without doubt one of the most explosive developments ever to have taken place in the telecommunications industry. At the same time, over the last years, we have seen the power of microprocessors double about every 18 months, becoming considerably smaller, cheaper and abundant; indeed, they are ubiquitous and are even finding their way into everyday objects. Those technology trendstiny, cheap processors with integrated sensors and wireless communications formed the technological basis for a new era of vast number of smart objects that communicate via wireless links.

Those incredibly diverse and complex wireless communication networks brought several fundamental technical issues: networks must be designed to carry out the intended functions in an efficient and predictable way; they must be manageable and upgradeable; and the most of all: reliable. Today, one of the major emerging threat against reliability is malware, that is, malicious self-replicating code. Threats posed by malware ranges from attacks against the confidentiality of the communication to attacks that actually alter the information traffic, hence destroying the integrity of the network [START_REF] Khouzani | Maximum damage malware attack in mobile wireless networks[END_REF].

Recently, malware outbreaks designed for personal computer environments on wired networks like those of Slammer and Code Red worms over the Internet have already inflicted severe economic damages, infecting thousands of hosts in short periods of time. Then, it is of vital importance to predict the limits of the damages that the attackers can inflict in large wireless networks, subject to its fundamental limitations, such as limited energy, unreliable communications and topology changes due mobility.

In this work we illustrate how the mean field approach can be used to reduce the complexity in the analysis of the damage that can be inflicted in a large network with opportunistic interaction of the objects. Opportunistic networking is one of the emerging communication paradigms in wireless mobile communications, in which the communication opportunities are based on sporadic and intermittent contacts. In contrast with the extensive work of malware propagation in wired networks ( [START_REF] Zou | The monitoring and early detection of internet worms[END_REF][START_REF] Garetto | Modeling malware spreading dynamics[END_REF][START_REF] Ganesh | The effect of network topology on the spread of epidemics[END_REF]), large wireless networks have not yet received similar attention, with a few exceptions ( [START_REF] Karyotis | Malware-propagative mobile ad hoc networks: Asymptotic behavior analysis[END_REF]6]). they give new insights to uncontrolled mean field framework which may be constrained (energy limitation). This helps in controlling the proportion of infected nodes.

Organization

The rest of the paper is structured as follows. In next section we present basics on stochastic games with random set of interacting players. In section 3, we overview existing mean field models in discrete time. We develop a controlled mean field framework in section 4 and in section 5 we provide a general convergence to mean field which is characterized by a stochastic differential equation and the payoff evolution are solution of partial differential equations. Finally, we apply the noisy mean field framework to opportunistic wireless large networks.

Basics of stochastic games with individual states

In this section we introduce basic notions of stochastic games and in subsequent sections we establish their convergence to the noisy mean field when the number of players goes to infinity under suitable conditions. A stochastic game with individual states and random set of interacting players is a collection Γ = N , S, (X j , Ãj , A j , r j ) j∈N , q, B n where • N is the set of players. The cardinality of N is n.

• S is a set of environment states. A state vector.

• For every player j, X j is its state space. The state of player j at time t is denoted by X n j (t) = (θ j , Y n j (t)) where θ j is the type, which is constant through the game. The set of possible states X j = {1, 2, ..., Θ} × Y j is finite. Y j may include other parameters, such as, space location, current direction and so on.

• For every player j ∈ N , Ãj is the set of actions of that player. A j : S × X j -→ 2 Ãj is a set-valued map (correspondence) that assigns to each state (s, x j ) ∈ S × X j the set of actions A j (s, x) that are available to player j. We denote the state-action space by SXA= (s, x, a):(s, x 1 , . . . , x n ) ∈ S× j∈N X j , a=(a j ) j∈N , a j ∈A j (s, x), ∀j ∈ N • For every player j ∈ N , r j : SXA -→ R is an instant payoff function for player j.

• q : SXA -→ ∆(S × j X j ) is a transition function where ∆(S × j X j ) is the space of probability distributions over S × j X j . The marginal of q is denoted by q j , which represents the transition probabilities of state of player j, given the state of all other players and the resource state.

• B n is a random set of players: B n ⊆ N denotes the set of interacting players at current time. The random selection will be explained in detail in the next section. The payoff function will be denoted by r B n j . How does the dynamic game evolves? The dynamic game starts at an initial (s 0 , x 0 ) ∈ S × j X j and is played as follows. Time space is N. A set B n t of players are randomly selected without replacement for a one-shot interaction. The current state s t , x t , B n t is known by all the selected players. At each time t ∈ N, each player j ∈ B n t , chooses an action a j,t ∈ A j (s t , x t ), and receives an instant payoff r j (s t , x t , a t ), where a t = (a j,t ) j∈N , and the game moves to a new state according to a probability distribution given by q(. | s t , x t , a t ) ∈ ∆(S × j ′ X j ′ ). The payoff and transitions of non-selected players are zero and the payoff function r B n t j of a selected player j does not depend on the state-actions of the non-selected players.

Strategies

Under complete information of past play, a history of length t + 1 corresponds to the sequence h t = (s 0 , x 0 , B n 0 , a 0 , s 1 , x 1 , B n 1 , a 1 , . . . , s t , x t , B n t ). A pure strategy σ j of player j is a mapping that assigns to every finite history h t an element of A j (s t , x t ). A mixed strategy is a probability distribution of the pure strategies. The set of mixed strategies is denoted by Σ j . A strategy profile (σ j ) j∈N is stationary if ∀ j, σ j (h t ) depends only on the current pair (s t , x t ). A stationary strategy of player j can be identified with element of the product-space of (s,x)∈S×X ∆(A j (s, x)). Every profile σ = (σ j ) j∈N ∈ j Σ j of mixed strategies, together with the initial state s 0 , x 0 induces a probability distribution P s 0 ,x 0 ,σ over the space of infinite plays SXA ∞ . We denote the corresponding expectation operator by E s 0 ,x 0 ,σ .

Payoff

We examine the finite-horizon payoff. Let σ be a mixed strategy profile. In the finite-horizon payoff, a player considers the cumulative payoff during the first T times. For every finite T ∈ N, the finite-horizon payoff for player j is

F j,T (s 0 , x 0 , σ) = E s 0 ,x 0 ,σ T -1 t=0 r B n t j (s t , x t , a t )1l {j∈B n t } (1)
where 1l {•} equals 1 if • is true and 0 if not.

The classic mean field approach

In this section we present a summary of the previous work on mean field and its relation with stochastic games.

Mean Field Interaction model [7]

The main interests of [START_REF] Benaïm | A class of mean field interaction models for computer and communication systems[END_REF] are models of interacting objects in discrete time, with finite number of states. The objects share local resources, which have a finite number of states. The objects are observable only through their state. In the limit, when the number of objects goes to infinity, it is found that the system can be approximated by a deterministic, usually non linear, dynamical system called the mean field limit. The mean field is in discrete or continuous time, depending on how the model scales with the number of objects. If the expected number of transitions per object per time slot vanishes when n grows, then the limit is in continuous time. Otherwise, the limit is in discrete time. One of the main advantages of the mean field approach is that it could drastically reduce the number of parameters in the analysis of large systems by using an aggregate dynamical system description.

X n (t) is called a mean field interaction model with n objects if: 1) X n (t) ∈ X n is an homogeneous Markov chain.

2) The transition kernel L n of X n (t) is invariant under any permutation of the labeling of the objects. That is,

L n (x, s; x ′ , s ′ ) = L n (x 1 , .., x n , s; x ′ 1 , .., x ′ n , s ′ )=L n (x p(1) , .., x p(n) , s; x ′ p(1) , .., x ′ p(n) , s ′ ) =P(X n 1 (t+1)=x ′ 1 , .., X n n (t+1)=x ′ n , S(t+1)=s ′ X n 1 (t)=x 1 , .., X n n (t)=x n , S(t)=s)
for any permutation p of the index set {1, 2, .., n}.

The occupancy measure is defined as the vector of frequencies of states:

M n x (t) = 1 n n j=1 1l {X n j (t)=x}
By the invariance assumption, (M n (t), S(t)) is an homogeneous Markov chain. Its state space is ∆ × S, where ∆ is the simplex {m = (m i ) :

i m i = 1, m i ≥ 0, i = 1,
∈ S, ∀x, x ′ ∈ X m ′ 1l { m-m ′ ≤2} (m x -m ′ x )(m x ′ -m ′ x ′ ) 1 δ n L n (dm ′ ; m, s) → 0 where L n (m ′ ; m, s) = P(M n (t+1)=m ′ |M n (t)=m, S(t)=s) H4: L n s;s ′ (m) is a smooth function of δ n and m. H5: f n (m, s) is a smooth function of δ n and m.
Later we will see that these assumptions are indeed relaxed in our approach. This theorem can be used to approximate M n (t) by the solution m(δ n t) of an ODE with the same initial condition, that is,

Theorem 1. If M n (0) → m 0 in probability [resp. in mean square] as n → ∞ then sup 0≤τ ≤T M n (τ ) -m(τ ) → 0 in probability [resp. in mean square], τ ∈ [tδ n , (t + 1)δ n ], where m(τ ) satisfies ṁ = f (m) m(0) = m 0
where f (m) = s w s (m)f (m, s) and w(m) is the invariant probability of the transition matrix L(m).

Mean Field Asymptotic of Markov Decision Evolutionary Games and

Teams model [START_REF] Tembine | Mean field asymptotic of markov decision evolutionary games and teams[END_REF] The main interests of [START_REF] Tembine | Mean field asymptotic of markov decision evolutionary games and teams[END_REF] are large populations of players in which frequent interactions occur between small numbers of chosen individuals. Each interaction in which a player is involved can be described as one stage of a dynamic game. The state and actions of the players at each stage determine an immediate payoff for each player as well as the transition probabilities of a controlled Markov chain. Each player wishes to maximize or minimize its expected payoff averaged over time. This model with a finite number of players, is in general difficult to analyze because of the huge state space required to describe all of the players. Then taking the asymptotics as the number of players grows to infinity, the whole behavior of the population is replaced by a deterministic limit that represents the system state, which is fraction of the population at each individual state that use a given action.

For large n, under assumptions analogous to H1-H5, the mean field converges to a deterministic measure that satisfies a non-linear ordinary differential equation for under any stationary strategy. They show that the mean field interaction is asymptotically equivalent to a Markov decision evolutionary game. The mean field asymptotic calculations for large n for given choices of strategies allows to compute the equilibrium of the game in the asymptotic regime.

Theorem 2. If M n (0) → m 0 in probability as n → ∞ then, for any stationary strategy u, and any time t, M n (t) converges in law to the solution of ṁ = f (u, m) m(0) = m 0 u is formally defined in the next section.

Controlled mean field interaction model

In this section, we introduce a controlled mean field interaction model. The finite version of this model is a particular case of stochastic games with individual states. We restrict our attention into a particular class of behavioral strategies within we are able to establish the mean field convergence. This restriction is due to the fact that when the number goes to infinity, the dimension of the set of stationary strategies goes to infinity as well as.

By letting the size of the system go to infinity, the discrete stochastic game problem is replaced by a limit of a system of Hamilton-Jacobi-Bellman equations coupled with a mean field limit ODE or coupled system of Bellman-Shapley optimality and discrete mean field evolution, that are deterministic and where the dimensionality of the original system has been transformed in the mass-behavior of the system. Time t ∈ N is discrete. The global state of the system at time t is (S(t), X n (t)) = (S(t), X n 1 (t), ..., X n n (t)). Denote by A n (t) = (A n 1 (t), . . . , A n n (t)) the action profile at time t. The system (S(t), X n (t)) is Markovian once the action profile A n (t) are drawn under Markovian strategies. We denote the set of Markovian strategies by U .

M n x (t) is the fraction of players who belong to the population of individual state x, as defined in section 3. Similarly, we associate the process

U n a (t) = 1 n n j=1 1l {A n j (t)=a}
to the fraction of actions.

Strategies and random set of interacting players

At time slot t, an ordered list B n t , of players in {1, 2, . . . , n}, is randomly selected without repetition as follows: First we draw a random number of players k t such that P(|B

n t | = k | M n (t) = m) =: J n k (m)
where the distribution J n k (m) is given for any n, m ∈ {0, 1 n , 2 n , . . . , 1} |X | . Second, we set B n t to an ordered list of k t players drawn uniformly at random among the n(n -1)...(nk t + 1) possible ones.

Each player such that j ∈ B n t takes part in a one-shot interaction at time t, as follows. First, each selected player j ∈ B n t chooses an action a j,t ∈ A(s, x j ) with probability u(a j | s, x j ) where (s, x j ) is the current player state. The stochastic array u is the strategy profile of the population.

Denote the current set of interacting players B n t = {j 1 , . . . , j k }. Given the actions a j 1 , ..., a j k drawn by the k players, we draw a new set of individual states (x ′ j 1 , ..., x ′ j k ) and resource state s ′ with probability L n s;s ′ (k, m, a), where a is the vector of the selected actions by the interacting players.

We assume that for any given Markovian strategy, the transition kernel L n is invariant by any permutation of the index of the players within the same type. This implies in particular that the players are only distinguishable through their individual state. Moreover, this means that the process M n (t) is also Markovian once the sequence of strategy is given. Denote by w n s,s ′ (u, m) the marginal transition probability between the resource states. Given any Markov strategy and any vector m of ∆, the resource state generates an irreducible Markov decision process with limiting invariant measure w s (u, m). Then, we can simplify the analysis by fixing the resource state S(t) = s without losing generality. The model is entirely determined by the probability distributions J n , the transition kernels L n and the strategy profile u.

Noisy mean field approach

We provide a general convergence result of the mean field to a stochastic differential equation and a martingale problem is formulated for the the law of the process M n t . We are able to establish a mean field convergence to non-deterministic differential equations, thus, extending the previous works in mean field interaction [START_REF] Benaïm | A class of mean field interaction models for computer and communication systems[END_REF][START_REF] Chaintreau | The age of gossip: Spatial mean-field regime[END_REF], in mean field Markov decision teams [START_REF] Tembine | Mean field asymptotic of markov decision evolutionary games and teams[END_REF][START_REF] Gast | Mean field for markov decision processes: from discrete to continuous optimization[END_REF] or in mean field Markov games [START_REF] Tembine | Mean field asymptotic of markov decision evolutionary games and teams[END_REF][START_REF] Weintraub | Oblivious equilibrium: A mean field approximation for large-scale dynamic games[END_REF]. We show that even if the expected number of players that do a transition in one time slot is not bounded, one can have a mean field limit, in such a case a stochastic one. This mean field limit is referred as noisy mean field.

Before presenting the main theoretical results of this paper, we first introduce some preliminary notions. The evolution of the system depends on the decision of the interacting players. Given a history h t = (S(0), X n (0), A n (0), . . . , S(t) = s, X n (t), A n (t)).

The random vector X n (t + 1) evolves according to the transition probability

L n (x ′ ; x, u, s) = P (X n (t + 1) = x ′ | h t ) . (2)
The term L n (x ′ ; x, u, s) is the transition kernel on X n under the strategy

U n . Let x n = (x n 1 , . . . , x n n ) such that 1 n n j=1 δ x n j = m and define L n (m ′ ;m, u, s) := (x ′ 1 ,...,x ′ n ) 1 n n j=1 δ x ′ j =m ′ L n (x ′ ; x, u, s) = P(M n (t+1)=m ′ |M n (t)=m, U n (t)=u, S(t)=s) = P(M n (t+1)=m ′ | ht )
where ht = (S(t ′ ), X n (t ′ ), A n (t ′ ), t ′ ≤t, S(t)=s, X n (t)=x n ).The term L n (m ′ ; m, u, s) corresponds to the projected kernel of L n . Below we give sufficient conditions on the transition kernels to get a weak convergence of the process M n t under the strategy U n (t).

We now present the main assumptions of this paper. A1: w s (u, m) is continuously differentiable in m and u. Note that this assumption was already implicit in the approach of [START_REF] Benaïm | A class of mean field interaction models for computer and communication systems[END_REF] via the smoothness of the drift. A2: There exists a continuous mapping f : R |X| × U × S -→ R |X| such that ∀s ∈ S;

lim n sup u∈U sup m ≤1 f n (m, u, s) δ n -f (m, u, s) = 0 where f n x (m, u, s)= m ′ ∈R |X| 1l { m ′ -m ≤2} (m ′ x -m x )L n (dm ′ ; m, u, s),
for x ∈ X , and s ∈ S. This assumption is analogous to H2 in which control parameters are added. A3: There exists δ n ց 0 and a continuous mapping a : R |X| × U × S -→ R |X|×|X| such that ∀s ∈ S;

lim n sup u∈U sup m ≤1
a n (m, u, s) δ n a(m, u, s) = 0 where

a n x,x ′ ,s (m, u, s) = m ′ ∈R |X| 1l { m ′ -m ≤2} (m ′ x -m x )(m ′ x ′ -m x ′ )L n (dm ′ ; m, u, s), (3)
for (x, x ′ , s) ∈ X 2 × S, where the third moment is finite. Note that under H3, a ≡ 0. A4: For all ǫ > 0 and ∀s ∈ S;

lim n sup u∈U 1 δ n m ′ ∈R |X| 1l { m ′ -m >ǫ} L n (dm ′ ; m, u, s) = 0
Note that assumption H3 and lemma 3 from [START_REF] Benaïm | A class of mean field interaction models for computer and communication systems[END_REF], implies A4. A4': ∀s ∈ S;

sup u∈U sup m∈R |X| sup n≥1 a n (m, u, s) δ n + f n (m, u, s) δ n < ∞
The smoothness assumptions H4-H5 in [START_REF] Benaïm | A class of mean field interaction models for computer and communication systems[END_REF] imply A4 ′ which implies A4 too.

Theorem 3. Let M n 0 -→ π in law where π is a probability measure. Under A1-A4 [resp. A1-A4'], the process M n t converges in law to a weak [resp. strong] solution of the Itô stochastic differential equation given by ,u,s) and w(m, u) is the invariant probability of the transition matrix L(m, u)

dm = f (m, u)dt + σ(m, u)dW t where σσ t = ã, f (m, u)= s w s (m, u)f (m, u, s), ã(m, u)= s w s (m, u)a(m
Proof. See Appendix.
Remark This result generalizes the deterministic mean field limit conditions established in [START_REF] Benaïm | A class of mean field interaction models for computer and communication systems[END_REF] for U equal to a singleton. It generalizes also the deterministic controlled mean field dynamics obtained in [START_REF] Tembine | Mean field asymptotic of markov decision evolutionary games and teams[END_REF] for stationary strategies.

Remark The conditions A1 -A4, A4 ′ are weaker than those given in [START_REF] Benaïm | A class of mean field interaction models for computer and communication systems[END_REF]. Under the conditions in [START_REF] Benaïm | A class of mean field interaction models for computer and communication systems[END_REF], the noise term 1 δn a n -→ 0 when n goes the infinity. Moreover, the continuity assumption on the drift limit is not needed. If f admits a unique integral curve, and ã is bounded and continuous then, the result applies as well. This allows us to apply it in wide range of networking scenarios with discontinuous drift limit but lower semi-continuity properties. Note that the uniformity in u may not be satisfied. In that case, a local mean field solution is derived. An example of such discontinuity is provided in [START_REF] Tembine | Mean field asymptotic of markov decision evolutionary games and teams[END_REF]. Our result extends also the convergence theorem in [START_REF] Tembine | Mean field asymptotic of markov decision evolutionary games and teams[END_REF] which was restricted to stationary strategies. Here u t is an admissible strategy at time t.

Definition 1 (Individual optimization framework). To a game as defined above, we associate a macroscopic population game, defined as follows. Each member j of the population, with state X j (t) and a population profile m[u](t). The initial condition of the game is S(0) = s, X j (0) = x, m[u](0) = m 0 . The population profile is solution to the noisy mean field evolution and S(t), X j (t) evolves as a jump process given by the marginal of q which depends on m(t), u(t) and the strategy u ′ (t) of j. Further, let V T (s, x, u ′ , u, m) be the T -stage payoff of player j in this game, given that S(0) = s, X j (0) = x and m(0

) = m 0 , i.e. V T (s, x, u ′ , u, m) = E g(m(T )) + T 0 r(S(t), X j (t), u ′ (t), u(t), m(t))dt s, x, m 0
where m is solution of the Itô SDE

dm(t) = f (m(t), u(t))dt + σ(m(t), u(t))dW t (4)
Using Itô's formula [START_REF]Stochastic Differential Equations: An Introduction with Applications[END_REF] the expected payoff evolution for a fixed horizon T is then given by

               v u ′ ,u (T, m) = g(m) -∂ t v u ′ ,u (t, m) = r(t, m, u, u ′ ) + x∈X fx (m, u)∂ mx v u ′ ,u (t, m) + 1 2 (x,x ′ )∈X 2 ãx,x ′ (m, u)∂ 2 mxm x ′ v u ′ ,u (t, m)
where v u ′ ,u (t, m) = E S,X j V T (s, x, u ′ , u, m) starting at time t and r(t, m, u, u ′ ) = E S,X j r(S, X j , u ′ , u, m) is the global expected mean field payoff.

Centralized mean field control: We now provide the feedback optimality principle for the global expected mean field payoff r(t, m, u, u). Please note that in this case, each member of the population uses the same strategy. A strategy u * (t) = φ(t, m) constitutes an optimal mean field solution if there exist a continuously differentiable function v(t, m) satisfying the following Hamilton-Jacobi-Bellman equation combined with the mean field SDE. By combining the system, one gets the mean field optimality: Proposition 1. The mean field optimality for horizon T is given by

                             v(T, m) = g(m) -∂ t v = sup ut∈Ut r(t, m t , u t , u t )+ x∈X fx (m t , u t )∂ mx v + 1 2 (x,x ′ )∈X 2 ãx,x ′ (m t , u t )∂ 2 mxm x ′ v   dm t = f (m t , u t )dt + σ(m t , u t )dW t , t > 0 m 0 = m where v = v(t, m).
Please note that this kind of systems may not have a solution. Counterexamples can be found in [START_REF] Tembine | Mean field stochastic games: Simulation, dynamics and applications[END_REF].

Discrete time mean field limit

If the step-size δ n is not vanishing when n goes to infinity, then the deterministic mean field limit is in discrete time and driven by the probability transition L t (u, m). Given an initial population profile m 0 and a terminal payoff, the sequence of population profiles {m t } t is driven by the transition probabilities {L t,s,x,s ′ ,x ′ (u, m t )} t .

m t+1 (x) = x ′ ∈X m t (x ′ )L t,x,x ′ (u t , m t ), (5)
where

L t,x,x ′ (u, m) = k≥0 s w s (u, m)L t,s,x,x ′ (u, m; k)J k (m), L t,s,x,x ′ (u, m; k)
is the limiting probability transition from x to x ′ when the resource state is s and the number of interacting players is k. Combining with the Bellman-Shapley optimality criterion, one gets the following system in the finite horizon case:

           v j,t (s, x, m) = max u j r(s, x, u, m t ) + s ′ ,x ′ L(s ′ , x ′ |s, x, u, m t )v j,t+1 (s ′ , x ′ ) m t+1 (x) = x ′ ∈X m t (x ′ )L t,x ′ ,x (u, m t ) (6)

Application to malware propagation

In this section, we apply the mean field approach to a controlled malware propagation in opportunistic networks. The malware propagation model is based on [START_REF] Benaïm | A class of mean field interaction models for computer and communication systems[END_REF] in which the impact of the control parameters is not examined and the player types are not used. The types can represent, for example, different operating systems, different versions of the operating systems, or patched/unpatched version of the same operating system. In most of the related work about malware spreading in large networks authors do not model the heterogeneity of systems which forms the network, as far as we know. This leads to a limitation in the results obtained, because different types of systems could lead for example to slower rates of propagation.

In this example we have mobile nodes that can be infected by a malicious code. There are two infected states: passive and active. Non infected nodes are susceptible. Then, the set of possible state of a node is {P, A, S} (for passive, active and susceptible) and the set of possible types is {θ 1 , θ 2 }. The state of the system at time t is X n (t) = (P 1 (t), P 2 (t), A 1 (t), A 2 (t), S(t)), where P (t) + A(t) + S(t) = n, j P j (t) = P (t), j A j (t) = A(t), ∀t and n is the total number of mobiles in the system. In this example, there is no resource. The occupancy measure is M n θ (t) = (P θ (t)/n, A θ (t)/n, S(t)/n) = (P n θ (t), A n θ (t), S n (t)). At every time step we want to control the proportion of infected nodes, which is I n (t) := A n (t) + P n (t). There are two fundamental ways to get infected:

1. Caused by a system flaw. (e.g: an exploit that could allow arbitrary code execution). 2. Caused by human flaw. (e.g: the user is deceived and executes a dangerous piece of code).

We can model this system as a controlled mean field interaction model. The interaction is simulated using the following rules: 1) A passive node may become susceptible (inoculation) with probability δ P .

2) A passive node with type θ may opportunistically encounter another passive node of type θ ′ , and both become active. This occurs with probability proportional to the frequency of other passive nodes at time t. For type θ, the probability is λ(P n θ ′ (t) -1 n 1l {θ=θ ′ } ). Note that the passive node can decide to contact the other passive node or not, so there are two possible actions: {m, m} (for meet and not meet). Those events will be modeled as a Bernoulli random variable with success (meeting) probability δ m , which represents u(m|P, θ). Here we model the possibility of getting infected by a system flaw.

3) An active node may become susceptible (inoculation) with probability δ A . 4) An active node of type θ may become passive with probability β

P n θ (t) h θ +P n θ (t) at time t.
Here is assumed that, at high concentrations of passive nodes, each active node infects some maximum number of passive ones per time step. This reflects finite total bandwidth. The parameter 0 ≤ β ≤ 1 has the interpretation of the maximum infection rate. The parameter 0 ≤ h θ ≤ 1 is the passive node density at which the infection proceeds at half of its maximum rate. Here we model the possibility of getting infected by a system flaw. 5) A susceptible node may become active with probability δ S . 6) A susceptible node may become passive via two ways. First, δ Sm is the probability of getting infected by a human flaw. In this case, the susceptible node can "decide" to get deceived or not, so there are two possible actions: {o, ō}. The stationary strategy in this case will be modeled as a coin toss with probability δ e . Second, η(P n θ (t) + P n θ ′ (t)) models the probability of encountering a passive node. In this case, the passive node can decide to contact the susceptible node or not, and it is modeled analogously to the other two cases.

At every time step, one of the transitions is randomly selected and performed. The number of nodes that do a transition in one time slot is always 0, 1 or 2. In order to control the infected population, each transition has a certain payoff contribution which could be 0 if no infected node is inoculated, 1/n if there is a node which is inoculated and -1/n if one node is infected. In Table 1 are the transition probabilities, the contribution to M n (t+1)-M n (t), the set of actions, and the contribution to the total payoff. Transition proba.

M n θ (t + 1) -M n θ (t) Actions Payoff 1 P n θ (t)δ P (-1, 0, 1)/n singleton set -1/n 2 P n θ (t)δ 2 m λ(P n θ (t) -1 n ) (-2, 2, 0)/n {m, m} 0 3 A n θ (t)δ A (0, -1, 1)/n singleton set -1/n 4 A n θ (t)β P n θ (t) h θ +P n θ (t) (-1, 1, 0)/n singleton set 0 5 S n (t)δ S (0, 1, -1)/n singleton set 1/n 6 S n (t)(δ e δ Sm + δ m ηP n (t))
(1, 0, -1)/n {o, ō, m, m} 1/n The intensity, that is, the probability that one arbitrary object does a transition in one time slot is of the order of 1/n. The drift, that is, the expected change of M n in one time step, given the current state of the system is:

f n θ (m) = nE(M n θ (t + 1) -M n θ (t)|M n (t) = m) =   -p θ δ P -2p θ δ 2 m λ np θ -1       

Homogeneous system

We briefly mention the homogeneous mean field. The drift is obtained by computing the expected changes in one time slot:

f n (m) = E(M n (t + 1) - M n (t)|M n (t) = m) = 1 n   -pδ P -2pδ 2 m λ np-1 n -aβ p h+p + s(δ e δ Sm + δ m ηp) 2pδ 2 m λ np-1 n -aδ A + aβ p h+p + sδ S pδ P + aδ A -sδ S -s(δ e δ Sm + δ m ηp)  
where m = (p, a, s). Then the limit is

f (m) =   -pδ P -2p 2 δ 2 m λ -aβ p h+p + s(δ e δ Sm + δ m ηp) 2p 2 δ 2 m λ -aδ A + aβ p h+p + sδ S pδ P + aδ A -sδ S -s(δ e δ Sm + δ m ηp)  
In all the simulations, we kept this parameters unchanged: β = 10 -2 , δ A = 5 × 10 -3 and δ P = δ S = δ Sm = 10 -4 . On parameter h (the passive node density at which the infection of active nodes proceeds at half of it maximum value) depends the stability of the system. Here we set h = 10 -2 in order to obtain an unstable behaviour. Regarding the control parameters, δ m = 1, δ e = 1 means no control.

We investigate the evolution of the system in the following scenarios:

• Trajectory of one run of the simulation (figure 1),

• Mean trajectory of multiple simulations (figure 2),

• Mean field limit trajectory (figure 3),

• Trajectory of the payoff function,

• Controlled mean field limit (figure 4),

• Trajectory of the heterogeneous malware propagation (figure 6),

• Optimal control under the mean field limit (figure 7),

• Noisy mean field (stochastic path), (figure 8).

These configurations are analyzed with control and without control parameters. We observe that the time mean 1 T T 0 m(s)ds is asymptotically close to the stationary point inside the limit cycle.

Uncontrolled behaviour

In figures 1,2,3 we can see the simulation results, using the well-known algorithm for exact simulation of a discrete time Markov chain. The initial configuration is (0.2, 0, 0.8). An oscillating behaviour of the total reward can be seen.

Controlled behaviour

In order to illustrate the control parameters, suppose we want to keep the proportion of infected nodes below 0.9 for all times. One simple way to achieve this is to reduce the contact tendency of a passive node. Then, we set δ m = 0.075, which is the more relevant control parameter. The results can be seen in figures 4 and 5.

Heterogeneous system

We investigate numerically the behavior of the mean field limit for two types θ and θ ′ . In figure 6 we can see that it is possible to stabilize the homogeneous system using classes.

Optimal strategy for the homogeneous system

Since the payoff function is the same for all the players i.e. r j (•) = r(•) the discounted stochastic game with common payoff can be transformed in a team problem. Moreover, the set of actions is the same for all the players. In figure 7, we plot the optimal strategy obtained by solving the system in the finite horizon case (equation ( 6)). The existence of a dominant strategy can be observed in the plot.

Noisy mean field

In order to give some feeling on how the noisy mean field evolves with time, we show in figure 8 two different realizations, for n = 1000 and n = 2000. The variance of the noise is obtained from σn = 1 δn f nf which its norm sup is bounded by 2Dλ n , then we have d m = f ( m)dt + σn dW t . Note that the smooth version of the last equation can be ṁn = f n (m n ) δn . It is worth to mention that, in this case, the simulation algorithm is only exact when the numerical time step vanishes.

Effect of the system size

In this subsection we examine the effect of n of the mean field object M n . We represent the evolution (in time) for different value of n and same state-component (see figures 9).

Concluding remarks

We have studied mean field stochastic games and established a mean field convergence to Itô stochastic differential equations. Our convergence result opens new questions for discontinuous mapping (namely the drift and the variance). In that case one may expect to have a stochastic differential inclusion as in [START_REF] Benaim | Stochastic approximation and differential inclusions[END_REF] for which the existence of solution needs more attention. This leads to ask what will happen when the discontinuities are coming from the strategies as it is the case for the best response correspondence? We do not have answers to this question and postpone it to future work. In practical scenarios, when the second moment conditions are not satisfied, the conditions for mean field convergence given in [START_REF] Tembine | Mean field asymptotic of markov decision evolutionary games and teams[END_REF][START_REF] Benaïm | A class of mean field interaction models for computer and communication systems[END_REF][START_REF] Gast | Mean field for markov decision processes: from discrete to continuous optimization[END_REF] are not applicable. In these cases, the weak conditions that we provided here cover a wide range of networking scenarios. The cycling behavior of the fraction of infected is sensitive to the control parameters. It would be useful to understand this relationship better by taking into account the network topology. [6] Y. Xu and W. Wang, "The speed of information propagation in large wireless networks," IEEE Infocom, 2008.

where E is a measurable subset in R |X | . We construct an interpolated process M n (.) from the discrete time population profile M n t with the step-size δ n . The time space is then the union of span in the form [t n k , t n k+1 ), where t n k = kδ n , k ∈ N. The martingale problem can be analogous defined as before.

Next, we construct a law π n,m,u of the new process with each interval in the form [δ n t, δ n (t + 1)]. It is not difficult to see that this covers the positive real line.

Then a martingale problem associated to π n,m,u and the transition kernel L n can be formulated. To do this, we define the martingale problem for the interpolated process from M n t denoted by M n (.) as above. We derive the generator of the process M n (.) and compute 1 δn L n . We show that the generator has limit if the assumptions A1-A4 are met. Let φ be an infinitely differentiable and bounded function defined on R |X | . Let O 12 (m, m ′ ) be the sum of first and second order in the Taylor development of φ i.e

O 12 := x (m ′ x -m x )∂ mx φ+ 1 2 x,x ′ (m ′ x -m x )(m ′ x ′ -m x ′ )∂ 2 mxm x ′ φ
Then, φ(m ′ , u)φ(m, u) -O 12 (m, m ′ , u) = O( m ′ -m 3 ). This means that if ψ n denotes the operator with the coefficients a n δn and f n δn i.e ψ n φ(m, u) = Under the above assumptions, we derive the tightness properties of the martingale law. Under the above assumptions, the sequence of processes {π n } n≥0 is pre-compact and any limit point as n -→ +∞ solves the martingale problem for f (., .) and a(., .) starting from m 0 . Finally, it follows from the assumptions A4' and the regularity of the coefficients f, a and the fact that the third order in Taylor approximation is bounded, the martingale problem has a unique solution determined by these properties up to a null set. Then, the noisy mean field limit (solution of the stochastic differential equation given by the generator) is the unique solution of the martingale problem, and π n -→ π the law of the process limit m. This completes the proof.

  The T -stage mean field optimization problem writes sup u Eg(m(T )) + T 0 r(t, m(t), u(t), u(t)) dt subject to the stochastic differential equation (4).

Figure 1 :

 1 Figure 1: Left: Path of one simulation without control. Right: Total reward.

Figure 2 :

 2 Figure 2: Leftmost: Mean of simulation trajectories without control. Rightmost: Total reward.

Figure 3 :

 3 Figure 3: Left: Deterministic mean field without control. Right: Total reward.

Figure 4 :

 4 Figure 4: Leftmost: Deterministic mean field with controlled passive node contact tendency. Rightmost: Total reward.

Figure 5 :Figure 6 :

 56 Figure 5: Evolution of the limit point of the time mean, for δ m = 0, 0.025, .., 0.925, 1.

Figure 7 :

 7 Figure 7: Optimal strategy using backward induction for the homogeneous system.

Figure 8 :

 8 Figure 8: Left: Path of one simulation of the noisy mean field, for n = 1000, Right: n = 2000.

Figure 9 :

 9 Figure 9: Mean field vs Simulation: Impact of n.

  m ′ , u)φ(m, u))L n (dm ′ ; m, u, s)ψ n φ(m, u, s) is bounded by c 1l { m-m ′ ≤1} mm ′ 3 L n (dm ′ , m, u, s) + 1l { m-m ′ >1} | φ(m, u)φ(m ′ , u) | L n (dm ′ , m, u, s).The last term is clearly bounded by2 φ ∞ 1l { m-m ′ >1} L n (dm ′ , m, u, s)which goes to zero by assumption A4. It remains to show that the first term goes to 0.1l { m-m ′ ≤1} mm ′ 3 L n (dm ′ , m, u, s) ≤ ǫ 3 + 1l { m-m ′ >ǫ} L n (dm ′ , m, u, s) -→ 0.

  Intensity vanishes at a rate δ n , that is, exists δ n , such that lim n→∞ δ n = 0 and

	2, ..|X |}. Since we are interested in the asymptotic results
	when n grows to infinity, we need to set the proper time scaling. Assume
	there exists a vanishing δ n called intensity, which can be understood as the
	probability that one arbitrary object does a transition in one time slot. Also
	define the drift f n (m, s) = E(M n (t+1)-M n (t)|M n (t) = m, S(t) = s), which
	is the expected change in the occupancy measure in one time slot. The main
	scaling assumptions in this model are:
	H1: The resource does not scale with n, that is, exists
	lim n→∞	L n s;s ′ (m) = L s;s ′ (m), ∀s, s ′ ∈ S
	where L n s;s ′ (m) = P(S(t + 1) = s ′ |M n (t) = m, S(t) = s). Moreover, the
	transition matrix L(m) has an unique invariant probability distribution.
	H2: lim n→∞	f n (m, s) δ n	= f (m, s), ∀s ∈ S

H3: Second moment of the number of object transitions per time slot is bounded, which is equivalent to ∀s

Table 1 :

 1 Probabilities, effects, actions and payoffs, (θ, θ ′ ∈ {1, 2}).

Appendix A. Sketch of the proof

Here we provide a proof sketch for noisy mean field convergence in law. We start by some notations. Given a test function φ defined on

since u is deterministic, φ depends only on m in the first and second terms.

With the notion of martingale, the question of characterization of the mean field trajectory can be formulate as martingale problem. One can ask if the property that D φ (.) is a martingale for all test functions φ uniquely characterizes the mean field m(.), up to an initial m(0).

The martingale problem for L z is : the existence of probability measure π defined on R |X | such that π(m 0 = m) = 1, and D φ (.) is a martingale for any test function. Is there at most one such π for each m?

We now define the martingale problem for the interpolated process from M n t denoted by M n (.).

L n (dm ′ ; M n δnt ′ , u, s), ∀t ′ ≥ 0.


Araks Martirosyan

Loïc Paulevé

Clair Poignard

Mireille Régnier

Jean-Marc Steyaert

Laurent Schwartz

A Dynamical Model for the Transmembrane Potential Regulation by pH

Three physico-chemical cell parameters have been shown by the biologists to be of importance for understanding the cell behavior: the transmembrane potential (TMP) which modifies substantially the cell functional activity, the intracellular pH i whose dynamics influences a number of other ionic fluxes and, to a lesser extent, the cell volume. In this work we present a dynamical model for understanding how the pH i influences the TMP. In this model we consider the HCO - 3 /Cl -and N a + /H + exchangers as well as the bicarbonate/carbon dioxide buffering mechanism and also the ions conductance dependence to the intracellular pH i .

We show that the conductances of the ions are the most important factors that have an impact on the TMP evolution and that the intracellular pH i controls the membrane ionic conductance.

Therefore, our new dynamical model support the hypothesis that pH i regulates the TMP dynamics due its influence on membrane ionic conductances.

INTRODUCTION

Most chemical reactions in our body are carried out by special proteins called enzymes. Moreover, the intracellular hydrogen ion concentration [H + ] (often denoted by

pH i = -log[H + ]
) is known to control these enzyme's activities, which happens to be optimal in a narrow range of H + concentration. Hence, understanding the hydrogen ions dynamics is a key for understanding the cell behavior. Two other physical cell parameters have been shown to be of importance by biologists: the cell volume whose too large fluctuations can endanger cell survival [START_REF] Lang | [END_REF] and the transmembrane potential (TMP) which modifies substantially the cell functional activity and affects its cell cycle [2][3][4][5][6][7].

Analysis of the experimental data [8][9][10][11] suggests that the intracellular pH i and the transmembrane potential are strongly correlated.

The purpose of the current work is to explicit this relationship by means of a dynamical model, which describes the phenomena. We then proceed to simulations in order to give a precise mathematical understanding of the relationship.

Related work

The relationship between the TMP and the ion fluxes has been addressed in several articles. In their seminal paper [3], Hodgkin and Huxley initiated the quantitative study of ion fluxes through the membranes applied to the conduction and excitation of nerve cells. In [2], Endresen et al. derive equations for the ionic currents that flow through channels, exchangers and electrogenic pumps. They establish the overall energy conservation of the cell and illustrate their theory in a model of spontaneously active cells in the cardiac pacemaker. More recently Clair et al. [12] have generalized the previous studies: their model links transmembrane potential, ionic concentrations and cell volume; they show that volume stabilization occurs within minutes of changes in extracellular osmotic pressure and infer a relationship between transmembrane potential and cell volume.

Contribution Our model is a generalization of the model described in [12], although we do not take into consideration the cell volume. It includes the HCO - 3 /Cl -and N a + /H + exchangers, the bicarbonate/carbon dioxide buffering mechanism and the pH i dependence of the ion permeabilities -a new concept motivated from various experimental studies [9,13,14]. As is now familiar, the mathematical model is given by a set of ordinary differential equations (ODEs). Numerical simulations show that the conductances of the ions are the most important factors that have an impact on the TMP evolution. As it is described in [9], the intracellular pH i controls the membrane ionic conductance and therefore we formulate the following hypothesis: intracellular pH i controls the TMP dynamics due to its influence on the membrane ionic conductances.

II. THE DYNAMICAL MODEL

In order to obtain a simple model of the pH influence on TMP through the permeabilities of ions, we consider hereafter that the outer concentrations of species is maintained constant -we assume here they are controlled by the external environment of the cell; and that the cell volume is constant.

A. Ion Fluxes and Transmembrane Potential

The transmembrane potential (TMP) sums up the electrical potentials due to the cell membrane and is driven by the ionic currents across the membrane. We use standard equations to derive TMP and ion fluxes evolutions [2,3,12] .

Given a ionic species S with valency (sum of charges) z S , its equilibrium potential E S can be obtained using the Nernst equilibrium equation:

E S = R.T F.z S log( [S] o [S] i ) ( 1
)
where R is the universal gas constant, F the Faraday constant, and T is temperature in Kelvin (fixed to 310K in the remaining of this paper).

[S] i (resp. [S] o ) denotes the inner (resp.) outer concentration of species S.

The resting TMP, denoted by E m , can then be expressed as the weighted sum of equilibrium potential of ions I, weights being their respective permeability through the membrane.

E m = S∈I P S .E S S∈I P S (2)
where P S is the membrane permeability for ion species S. In our model I = {K + , N a + , Cl -}.

Given a TMP value V m , the ionic current I S for S is then determined by first deriving its conductance g S through the membrane (Eq. ( 3)), and is proportional to the difference between V m and the equilibrium potential of S (Eq. ( 4)). The evolution of the TMP is obtained with the Kirchoff law (Eq. ( 5)).

g S = P S .z 2 S . F 2 2.R.T . [S] o .[S] i (3)
I S = g S (V m -E S ) (4)
dV m dt = -S∈I I S C m (5)
where C m stands for the membrane conductance (fixed at 0.07 F.m 2 ).

Finally, the evolution of the inner concentration of species S due to its current follows the following equation:

d[S] Flux i dt = -I S A z S .F.V (6)
where A and V are respectively the surface area and the volume of the cell. In the rest of this paper, the volume is assumed constant and the ratio A/V has been fixed to 3/r, with r = 10 -5 m being the radius of the cell (assumed spheric).

B. pH-Dependent Permeabilities

The existence of a pH-dependence for ions permeabilities has been hypothesized in several experimental studies [9,13,14]. The proposed equation is inspired from [9], where it has been successfully fitted to experimental data.

In the scope of this article, we focus on having only Na + permeability being pH-dependent, hence time-dependent, with [H] t i being the H + concentration at time t:

P N a + (t) = P N a + 10 -S P .pK N a + 10 -pK N a + .S P + [H] t i S P (7)
where pK N a + is the midpoint of pH-dependence, that is

P N a + (t) = P N a + 2
when [H] t i = 10 -pK N a + . S P is the stiffness of this pH-dependence. Fig. 1 depicts the evolution of P N a + (t)/P N a + with different parameters. In such a setting, P N a + should be considered as the maximum Na + permeability.

The above equation can be straightforwardly generalised for introducing pH-dependence for each ion permeabilities.

As it is discussed in Subsect. III B, the result on TMP regulation by pH relies on the varying ratios between ions permeabilities, where the effect is modulated by inner and outer ionic concentrations. In the scope of this article, we selected parameters such that a varying (plain curve) pK N a + = 6.9, S P = 4; (dash-dotted curve) pK N a + = 6.9, S P = 2; (dotted curve) pK N a + = 6, S P = 2.

Na + permeability has the greatest impact on the TMP. With other concentrations and permeabilities settings our model should provide different conclusions.

C. pH Regulation Mechanisms

The regulation of pH within the cell is a very complex machinery involving numerous various biological processes. We consider here three important mechanisms: the bicarbonate buffer, and electro-neutrals Na + /H + and HCO - 3 /Cl -ionic exchangers. a. Carbon dioxide/bicarbonate buffer Although CO 2 does not contain hydrogen ions it rapidly reacts with water to form carbonic acid H 2 CO 3 , which further dissociates into hydrogen and bicarbonate ions HCO - 3 . For the sake of simply, we retain a single reversible reaction:

CO 2 + H 2 O k 1 -- k 2 HCO - 3 + H + (8)
where the ratio k 2 /k 1 is the acid dissociation constant pK a = 6.37. The water has to be considered as as abundant as needed.

Hence, this reversible reaction acts as a H + buffering mechanism, leading to a fixed balance between CO 2 and HCO - 3 •H + concentrations. The impact of this buffering on species concentrations evolution can be expressed as follows.

d[HCO 3 ] Buf i dt = d[H + ] Buf i dt = - d[CO 2 ] Buf i dt = -k 1 [HCO - 3 ] i [H + ] i + k 2 [CO 2 ] i (9)
b. Na + /H + unidirectional exchanger The role of Na + /H + exchanger for limiting the cell inner acidity has been underlined by several experimental studies [15][16][17][18][19] Basically, the exchanger is fully active when [H + ] i is above a certain threshold [H + ] * .

In addition, this exchanger can be stopped when the cell obtain a too large volume. As we consider a constant cell size in our model, we ignored this latest feature, and propose the following equation for H + evolution concentration, assuming external concentrations constant:

d[N a + ] nhe i dt = - d[H + ] nhe i dt = k N HE 1 + tanh(α([H + ] i -[H + ] * ) + tanh -1 (1 -)) 2 ( 10
)
where k N HE is the maximum rate of the exchanger and tanh the hyperbolic tangent. This term will tend toward k N HE when [H + ] i [H + ] * , and tend toward zero when [

H + ] i ≤ [H + ].
α is a rescaling factor for the H + concentration difference, and has been set to 10 10 in our simulations.

c. HCO - 3 /Cl -bidirectional exchanger Due to the bicarbonate buffering, inner pH can also be controlled by varying HCO - 3 concentration. According to [18,19], the HCO - 3 /Clm exchanger also play a role in pH regulation by acting on HCO - 3 inner concentration in order to maintain the pH around a certain threshold. For the sake of simplicity, we set this H + threshold identical to Na + /H + mechanism (while it appears to be lower in practice [19]).

Assuming external concentrations constant, the activity of this exchanger is limited by inner HCO - 3 and Cl -inner concentrations, and its direction depends on the sign [H + ] i -[H + ] * . When [H + ] i is above (resp. below) the threshold, the exchanger makes enter (resp. exit) HCO - 3 against Cl -. This can be summarized by the following equation:

d[HCO - 3 ] ae i dt = - d[Cl -] ae i dt = k AE .tanh(α([H + ] i -[H + ] * )) [HCO - 3 ] i β 1 + [HCO - 3 ] i [Cl -] i β 2 + [Cl -] i (11)
where k AE is the maximum rate of the exchanger, and β 1 and β 2 the respective saturation constants for HCO - 3 and Cl -(set resp. to 0.001 and 0.03).

D. Full Model

In addition of the dynamics described above, we also consider a (constant) metabolic production of H + and CO 2 . The full model can then be summarized as follows, where Na + permeability is computed using Eq. (7).

d[H + ] i dt = H prod + d[H + ] Buf i dt - d[H + ] nhe i dt d[CO 2 ] i dt = CO 2prod - d[CO 2 ] Buf i dt d[N a + ] i dt = d[N A + ] Flux i dt + d[N A + ] nhe i dt d[Cl -] i dt = d[Cl -] Flux i dt - d[Cl -] ae i dt d[K + ] i dt = d[K + ] Flux i dt d[HCO - 3 ] i dt = d[HCO - 3 ] Buf i dt + d[HCO - 3 ] ae i dt dV m dt = -S∈I I S C m III. RESULTS
We first show simulation results supporting the pH influence on the dynamics of TMP for particular permeabilities of K + , Na + and Cl -, and having permeability of Na + being pHdependent. Then, the sensitivity of this influence with respect to permeabilities parameters is studied in Subsect. III B. In particular, the simulations strongly suggest that the Na + pH-dependence plays a crucial role, and that the ratio of actual Na + permeability against K + permeability has to be above a certain threshold.

The selected parameters for our model are listed in Table I. Concentrations and permeabilities have been chosen to be within classical ranges [12,19,20]. For the purpose of our study, we enforce the pH to increase from 6.8 to 7.2 (i.e. [H + ] i is divided by ≈ 2.5); which corresponds to some observed pH variation amplitudes during cell cycles [4,5].

Simulations have been performed using the stiff ODE integrator shipped with Scilab [START_REF]Scilab Enterprises, Scilab: free open source software for numerical computing[END_REF].

A. pH influences the dynamics of TMP

In order to appreciate the effect of the pH-dependence of ions permeabilities on TMP dynamics, we first did some experiments without this pH-dependence enabled.

Fig. 2 shows the integration of our model including the ion fluxes and TMP equations (Subsect. II A) for the dashed curve; plus the pH regulation mechanisms (Subsect. II C) for the dash-dotted curve. In both cases, the permeability of Na + has been fixed to ≈ 0.28.P N a + , which corresponds to its pH-dependent value when pH is 6.8 (Subsect. II B).

We observe that the regulation of pH has no impact on TMP dynamics. Indeed, despite the activation of Na + /H + and HCO - 3 /Cl -exchangers, the concentration of H + being very small in front of Na + and Cl -, the variation induced in Na + and Cl -concentrations is actually negligible. Thus, the sole regulation of pH should not impact the TMP.

We also remark that the TMP is not stable in our model settings. This is expected regarding the resting potential of K + and Na + ions (E K + ≈ -95mV and E N a + ≈ +68mV) their gradients make them inherently decrease and increase, respectively. As discussed in Sect. IV, future work may develop a more complete model, including notably the Na + /K + .

6.8

6.9 The integration of our complete model compared to the sub-model having only TMP equations is presented in Fig. 3 and Fig. 4.

The pH-dependence of Na + permeability makes it to modify the TMP significantly when pH varies. In the first 40 seconds, while pH moves from 6.8 to 7.2, the permeability of Na + goes from ≈ 0.28.P N a + to ≈ 0.94.P N a + . Hence, Na + current is strongly increased, impacting the global TMP. After 40 seconds, pH is stabilized, Na + permeability becomes constant and remain at a high value compared to the initial one, inducing a greater constant TMP decrease.

Fig. 4 details the concentration evolution of K + , Cl -, Na + , and HCO - 3 in same setup as in Fig. 3. As previously discussed, the large variation of K + and Na + dynamics should not be attributed to the activation of Na + /H + exchanger, but to the evolution of TMP (for K + ) and permeabilities (for Na + ). The tiny variation of Cl -and HCO - 3 dynamics is mainly due to H + regulation by Cl -/HCO - 3 exchanger.

B. Parameters sensitivity

As recalled in Subsect. II A, at constant temperature, while the equilibrium potential of ions is solely determined by the ratio of their inner and outer concentration, the resting TMP follows a weighted sum of those equilibrium potential, the weights being the ions permeabilities through the membrane. The flux of each ion is then proportional to the difference between its proper equilibrium potential and the global TMP, and proportional to its permeability.

In such a setting, the TMP is controlled by both ions inner and outer concentrations and their permeabilities. Focusing on permeabilities, the crucial factor is their mutual ratios (because of the weighted sum). Hence, the pH influence on TMP is explained by the ratio between Na + and both K + and Cl -permeabilities which significantly varies alongside the H + concentration.

Given the selected initial ions concentrations (Table I), the resting TMP is around -91mV, and the ions resting potential are E K + ≈ -95mV, E N a + ≈ +68mV and E Cl -≈ -94mV.

Hence, we notice that E K + and E Cl -are actually close to initial TMP, resulting in low ionic currents through the membrane. On the contrary, E N a + present an opposite value, indicating a high current, however restrained by the (typically) low Na + permeability. But, because E N a + is so different from E K + and E Cl -, even with a low permeability with respect to K + , a variation in Na + permeability should have impact a noticeable impact on the TMP. This is notably illustrated in Fig. 5 by integrating TMP dynamics with different permeability ratios settings (left figure), or with different pH-dependence for Na + permeabilities.

The other parameters are identical as in previous section, and the full model is considered. (right) Integrated TMP with S P = 0.1 (dash-dotted), 2 (dashed), 4 (plain).

-95m V -90m V -85m V -80m V -75m V 0 20 40
In Fig. 5(left), the maximal permeability P N a + ranges through 0.004P K + , 0.04P K + (default value), and 0.07P K + , showing the strong impact of the P N a + /P K + ratio. We notice that a too small ratio (below 1%) makes E N a + negligible for TMP.

In Fig. 5(right), the stiffness of the pH-dependence for P N a + is controlled by tuning the S P parameter. With a low S P parameter, the Na + slowly varies with respect to pH, while a high S P induces a large variance. As predicted, the higher S P , the stronger TMP regulation.

IV. DISCUSSION

We presented a new dynamical model of ions fluxes through the cell membrane which emphasizes the regulation of the cell transmembrane potential (TMP) by the varying intracellular pH and ions permeabilities.

As stated in the introduction and in Subsect. II B, various experimental studies already support the hyothesis that inner pH controls ions permeabilities, and thus affect the cell TMP. Our model brings a mathematical analysis of this possibly involved regulation mechanism, and calls for further experimental validations.

Future work will incorporate additional ions regulation mechanisms, such as the Na + /K + pumps [START_REF] Armstrong | [END_REF], and will consider a dynamical cell volume with osmolarity pressure (as it has been done in prior work in [12]). Finally, the modelling of the impact of pH on

Figure 1 :

 1 Figure 1: Proportion of the maximal permeability P N a + in function of pH. (plain curve) pK N a + =

Figure 2 :

 2 Figure 2: pH (left) and TMP (right) dynamics with only TMP equations (dashed curve) and with TMP together with pH regulation mechnisms (dash-dotted curve).

Figure 3 :Figure 4 :

 34 Figure 3: pH (left) and TMP (right) dynamics with only TMP equations (dashed curve) and with the full proposed model (plain curve).

Figure 5 :

 5 Figure 5: (left) Integrated TMP with P N a + = 0.004 (dashed), 0.04 (plain), 0.07 (dash-dotted);


Sébastien Guigard

Nathalie Henrich Bernardoni

Jocelyne Sarfati

Christophe Savariaux

Philippe Chaffanjon

Assessment of vocal morbidity following cervical endocrine surgery