Arran Hodgkinson
email: arran.hodgkinson@googlemail.com

Ovidiu Radulescu
email: ovidiu.radulescu@umontpellier.fr

Gilles Uzé

Dumitru Trucu
email: trucu@maths.dundee.ac.uk

Giles Uzé
email: giles.uze@umontpellier.fr

Giles Uzé

Signal Propagation in Sensing and Reciprocating Cellular Systems with Spatial and Structural Heterogeneity

Keywords: Mathematics Subject Classification (2000) AMS Subject Classification: 22E46, 53C35, 57S20 Keyword1, keyword2, keyword3

HAL is

1 Introduction

Sensing and Reciprocating Systems and their Mathematical Treatment

In order for biological systems to initiate changes in behaviour at the scale of a group of cells or of a tissue in response to a localized event, it is necessary for small signals to be transformed into large signals and sequentially communicated to other cells. This is no more apparent than in the human immune response where T-cells are actively recruited to the site of infection through the amplification and dispersion of the precursor signal [START_REF] Mackay | Chemokine receptors and t cell chemotaxis[END_REF]. The intermediate signal must be received and amplified, in order that distant cells may receive the signal with sufficient veracity as to respond.

In the case of the immune system, the cell-to-cell communication can be at least partially orchestrated by dynamic changes of the cell membrane receptors and by secretion of communication proteins such as chemokines [START_REF] Mackay | Chemokine receptors and t cell chemotaxis[END_REF] and cytokines [START_REF] Oyler-Yaniv | A Tunable Diffusion-Consumption Mechanism of Cytokine Propagation Enables Plasticity in Cell-to-Cell Communication in the Immune System[END_REF]. Other cell-to-cell communication and amplification mechanisms are used by bacteria in a phenomenon known as "quorum sensing" [START_REF] Ng | Bacterial Quorum-Sensing Network Architectures[END_REF] and by yeast to optimize mating efficiency [START_REF] Barkai | Protease helps yeast find mating partners[END_REF]. In order to synchronise the phenotypes expressed by a local group of cells, bacteria and yeast posses internal feedback loops that amplify incoming diffusible chemical signals. Similar examples where local behavior spreads by cell-to-cell communication can be found in animal development, when blocks of tissues can be developed from sheets of cells by a phenomenon called "community effect" [START_REF] Gurdon | A community effect in animal development[END_REF] or when cell fate is specified by "sequential patterning" such as in the spatial regulation of Delta-Notch signalling [START_REF] Hoyos | Quantitative variation in autocrine signaling and pathway crosstalk in the caenorhabditis vulval network[END_REF][START_REF] Henrique | Maintenance of neuroepithelial progenitor cells by delta-notch signalling in the embryonic chick retina[END_REF]. Collective synchronous behaviour of cells is also needed in insulin secretion by pancreatic islets but, in this case, the possible cell-to-cell communication mechanisms are still under debate [START_REF] Pedersen | Intra-and inter-islet synchronization of metabolically driven insulin secretion[END_REF]. We call such systems sensing and reciprocating systems (SARs), on the basis that the initial chemical signals are replicated and amplified, which is similar to the concept of secrete and sensing cells [START_REF] Maire | Molecular-Level Tuning of Cellular Autonomy Controls the Collective Behaviors of Cell Populations Article Molecular-Level Tuning of Cellular Autonomy Controls the Collective Behaviors of Cell Populations[END_REF][START_REF] Olimpio | Out-of-equilibrium statistical dynamics of spatial pattern generating cellular automata Kavli Institute of Nanoscience[END_REF].

SARs are ubiquitous in biology and some mathematical models dealing with properties of such systems exist. The versatility of collective properties of secrete and sensing cells was studied using phenomenological, compartment based models and ordinary differential equations (ODEs) [START_REF] Youk | Secreting and sensing the same molecule allows cells to achieve versatile social behaviors[END_REF]. The same type of formalism was used for metabolic synchronisation of insulin secretion in islets [START_REF] Pedersen | Intra-and inter-islet synchronization of metabolically driven insulin secretion[END_REF] and for studying cell-to-cell communication in the immune system [START_REF] François | The Case for Absolute Ligand Discrimination : Modeling Information Processing and Decision by Immune T Cells[END_REF]. ODE based models allow rather detailed descriptions of intracellular signalling and metabolic dynamics but do not cope accurately with cell proliferation, migration, and cell-to-cell interactions.

Although not yet used for SARs, frameworks based on partial differential equations (PDEs) could integrate many of these processes and explain aspects related to spatial heterogeneity such as the role played of spatial arrangement of cells in determining the conveyance of these signals [START_REF] Oyler-Yaniv | A Tunable Diffusion-Consumption Mechanism of Cytokine Propagation Enables Plasticity in Cell-to-Cell Communication in the Immune System[END_REF]. However, in PDE models, non-spatial heterogeneity, resulting from the fact that cells in close spatial proximity do not necessarily respond synchronously to stimuli, is lost by averaging. This "structural heterogeneity" can be an essential part of a complex cell dynamics, in which cell sub-populations behave differently to the average, and may be essential to understanding the complicated dynamics of biological systems. As an example, such models have predicted that below a certain threshold value interferon signalling allows the activity of the cellular population to decay entirely [START_REF] Hart | Paradoxical signaling by a secreted molecule leads to homeostasis of cell levels[END_REF].

A paradigm which seems appropriate to exploring the possible structural dimensions of biological problems, in a mathematical context, is that of the continuous structural approach [START_REF] Bekkal Brikci | An age-and-cyclinstructured cell population model for healthy and tumoral tissues[END_REF][START_REF] Lorz | Populational adaptive evolution, chemotherapeutic resistance and multiple anti-cancer therapy[END_REF][START_REF] Lorz | Modeling the Effects of Space Structure and Combination Therapies on Phenotypic Heterogeneity and Drug Resistance in Solid Tumors[END_REF]. This approach encompasses the genetic or epigenetic state of a cell, under temporal conditions which are consistent with the continuous nature of dynamic biological problems by employing the application of PDEs in structure, rather than in spatial position. On the other hand, these approaches neglect the spatial dimensions associated with chemical communication between cells and, thusly, do not provide the descriptive breadth necessary to analyse these situations.

One recent "spatio-structural-temporal" (SST) framework, which demonstrates the potential to represent greater details of dynamical processes in dimensions of both structure and space, was developed in order to model the urokinase plasminogen activator system in breast cancer [START_REF] Domschke | Structured models of cell migration incorporating molecular binding processes[END_REF][START_REF] Trucu | Multiscale computational modelling and analysis of cancer invasion[END_REF]. Herein, we present a similar derivation in order to augment the generality of this framework and present a modelling form capable of capturing the intricacies, and important heterogeneous features of SARs. Compared to [START_REF] Domschke | Structured models of cell migration incorporating molecular binding processes[END_REF][START_REF] Trucu | Multiscale computational modelling and analysis of cancer invasion[END_REF] we introduce new metabolic structural variables and conjugated advection fluxes that are derived from the continuity equation and Liouville's theorem. These variables are needed for modelling stimulated amplification in SARs. The use of Liouville theorem is a major advance in the SST framework as it can relate any single cell ODE dynamics to population dynamics in structure space.

An example of SAR system: Cellular Interferon System

We look, here, specifically at a detailed model for the interferon (IFN) binding process of a given cell and the concurrent metabolic processes that result from this binding process. This SAR shall serve as an exemplar biological system on which to base models that will explore the efficacy of the framework to be proposed.

There are two primary forms of INF, of interest, namely INFα and INFβ. Their ability to activate a cell's internal infrastructure is dependent on their ability to concurrently bind the IFN-α/β receptors 1 (INFAR1) and 2 (IN-FAR2) on the surface of the cell. The association rate of IFN with IFNAR2 is approximately 10× that of IFNAR1, therefore the primary interaction is with the tyk2 signalling complex of IFNAR2 [START_REF] Gavutis | Lateral Ligand-Receptor Interactions on Membranes Probed by Simultaneous Fluorescence-Interference Detection[END_REF][START_REF] Gavutis | Determination of the Two-Dimensional Interaction Rate Constants of a Cytokine Receptor Complex[END_REF]. The essential IFN signalling protein, however, is the Jak1 tail of IFNAR1 and so IFNAR1 is preferentially recruited to the location of the bound IFN molecule [START_REF] Gavutis | Lateral Ligand-Receptor Interactions on Membranes Probed by Simultaneous Fluorescence-Interference Detection[END_REF][START_REF] Gavutis | Determination of the Two-Dimensional Interaction Rate Constants of a Cytokine Receptor Complex[END_REF]. These tyk2 and Jak1 protein phosphorylate one another to initiate what is known as the Jak-Stat pathway [START_REF] Stark | How cells respond to interferons[END_REF].

The Jak-Stat pathway is predicated on the fact that the phosphorylated Jak1-tyk2 complex is capable of phosphorylating the transcription factors Stat1 and Stat2. These two factors are then able to phosphorylate and bind the interferon regulatory (transcription) factor 9 (IRF9) in order to form the IFN stimulated gene factor 3 (ISGF3) complex [START_REF] Stark | How cells respond to interferons[END_REF][START_REF] Samuel | Antiviral Actions of Interferons.pdf[END_REF], which is capable of entering the nucleus [START_REF] Lin | Virus-Dependent Phosphorylation of the IRF-3 Transcription Factor Regulates Nuclear Translocation, Transactivation Potential, and Proteasome-Mediated Degradation[END_REF]. Once having achieved this step, this complex can bind to the promoter region of IFN stimulated genes (ISGs) and effectively initiating trancription of these ISGs [START_REF] Stark | How cells respond to interferons[END_REF][START_REF] Samuel | Antiviral Actions of Interferons.pdf[END_REF].

In turn, transcription of ISGs is responsible for the downstream phosphorylation of IRF7 and IRF3, who are directly responsible for the promotion of IFN-α and IFN-β genes [START_REF] Haller | The interferon response circuit : Induction and suppression by pathogenic viruses[END_REF]. Therefore, this system can be looked at through the simplified lens of two major and important processes:

(a) the binding of interferon to the surface of the cell, and (b) the activation of the metabolic pathway which eventually leads to the creation of new IFN molecules.

We use the phrase 'metabolic activation' in order to characterise the state of the cell in terms of the chemical activity levels of those proteins involved in the Jak/Stat pathway and, ultimately, the transcription of the genes necessary for the synthesis of IFN. Thus, when one describes the metabolic activation of the cell, with regards to the IFN pathway, one is actually describing, in some way, the spatially differentiated presence of IRF7/IRF3 within the cell (Figure 2). Moreover, one review of experimental data plotted the relationship over time between the activation of genes within the cell and the fractional levels of bound and unbound surface receptors, for both INFα2 and INFβ [START_REF] Schreiber | The molecular basis for functional plasticity in type I interferon signaling[END_REF]. This graph importantly showed that, for low levels of INFα2, as the number of surface receptors decreased, the metabolic activation level rose concurrently. Further, as genetic activation levels decreased, one could observe a corresponding normalisation of the fractional surface receptor levels [START_REF] Lavoie | Cytokine Binding and activity of all human alpha interferon subtypes[END_REF]. Comparably, for high levels of INFβ, one finds that the cells genetic mechanism is activated in a locally irreversible process and that the fraction of IFNAR1 receptors is maintained at approximately 40% [START_REF] Schreiber | The molecular basis for functional plasticity in type I interferon signaling[END_REF].

In order to demonstrate the descriptive power within the existing modelling frameworks, we choose the biological interferon system in T-cells as an illustrative example of such a system of SARs. This will serve as a comparative case for the development of a framework, which is capable of significantly improving upon one's existing capacity.

A Simple, Continuous Mathematical Model for the IFN System

If one were to create the simplest possible system of SARs, one would begin with only the population of SARs, themselves, and the molecular population of SAR diffusing ligands. In reality, however, these systems are rarely as simplistic as this and often require consideration of spatially intermediate cells which may mediate the levels of the SAR ligands, by consuming these proteins without reciprocally producing them. This is the case in the biological IFN system and, as such, we call such intermediate cells 'consumers' and the SAR cells as 'producers', within a system that considers only such a responsive protein.

Therefore, begin by defining a temporal domain, given by I = [0, T ] with t ∈ I, and a two-dimensional spatial domain, given by D ⊆ R 2 with x ∈ D. We then write cellular population functions such that c 1 : I × D → R gives the population of interferon producing cells and c 2 : I × D → R gives the global population of consumer cells, whilst m : I × D → R gives the nondimensionalised concentration of IFN molecules.

In order to write as simple a model as is possible, we begin by ignoring all dynamics in the cellular populations are given simply by c 1 (t, x) := c 1 (0, x) and c 2 (t, x) := 1, respectively. This is so that one might analyse only the communicative capabilities of the interferon itself.

We then write the dynamics of the system as a whole as an ordinary differential equation in m(t, x), such that the spatial dynamics are given entirely by the diffusion of this molecule in the solution. Interferon is then systematically consumed by c 2 , at a rate λ, and is autoreplicated within c 1 cells, at a rate φ 2 , and where this autoreplication is further stabilized by negative self-regulation, with the rate constant φ 3 . Therefore, we have that

dm dt = D m ∇ 2 x m -λmc 2 + (φ 2 m 2 -φ 3 m 3 )c 1 (1)
where D m is the coefficient for diffusion of interferon. Simulations were performed for this system using a 4 th order Runge-Kutta predictor and MacCormack corrector, with a central difference formula used for the calculation of diffusion terms. Initial conditions for the producer cells are given by

c 1 (0, x) = 5 j=1 5 i=1 exp -(x 1 - 1 2 i) 2 -(x 2 - 1 2 j) 2 ,
and for the interferon concentration is given by the Gaussian distribution

m 0 := m(0, x) = exp -(x 1 - 1 2 ) 2 -(x 2 - 1 2 ) 2 ,
with the rate constants given by D m = 10 -3 , φ 2 = 3 4 , and φ 3 = 1 8 and λ being variable between simulations.

The results for the simulation of system (1) show, most simply, that communicative capability increases with decreasing values for affinity of interferon for its consumer cells (Fig. 1). The approximate threshold value for which this is true falls in the interval λ ∈ (0.1, 0.15) 1, given the values chosen for D m , λ, φ 2 , φ 3 .

This may, to some extent, give a mathematical explanation for why it may be biologically advantageous to maximise the utilisation of lower affinity interferon in a system where one wishes to stop the spread of the infection. It could be that cells employ this methodology in order to spread a panic signal upon the initial detection of a virus and initialisation of a local interferon signal.

The explanation given by this simple model, however, does not explain to one the nature of the interaction between molecules and cells that allows this system to proffer communicative capabilities as it does. For example, one artificially introduces the notion that increasing the affinity of IFN molecules will increase their consumption but one is left wondering what effect this alteration should have on the interaction with producer cells. It is difficult to intuit, also, how this increase in affinity should change the interactions that impact the metabolism of IFN within the cell. One might expect that affinity would increase production but would it also increase feedback sufficiently to dampen that response? Alterations to equation 1, however, require suppositions on the desired final behaviour of the system, rather than a priori biological assumptions.

Therefore, in response to this fundamental issue, we aim to create a more biologically descriptive model that will serve to quantify dynamics in the cellsurface receptors; the binding of these receptors by free molecules; and the consequential alterations in metabolism in a spatial context. We will also realise the interactions between these various dynamical behaviours, in order that one might better understand how the biological reality is affected by changing individual characteristics.

General SAR model within the SST framework

We introduce here a general SST model for SAR systems. Various instances of this model can serve to study different biological problems.

In this framework we consider that cells of the same type can differ in their states. The cell state is described by three variables ξ ∈ Υ ⊂ R υ , y ∈ P ⊂ R p and α ∈ Γ ⊂ R γ , where ξ, y, α represent the total density of receptors on the cell membrane; the part of receptors that have bound ligands; and the metabolic variables, respectively. We consider that there are q different diffusible ligands of concentrations m l (t, x), 1 ≤ l ≤ q. As a simplifying assumption we consider that ligands m k bind with no competition to their cognate receptors ξ i , 1 ≤ i ≤ υ. Competition could be easily introduced by considering A spatially and structurally heterogenous cell population is described by a structured cell density, namely by a positive, integrable function ĉ(t, x, ξ, y, α), with t ∈ (0, T ], x ∈ D ⊂ R d , ξ ∈ Υ , y ∈ P, and α ∈ Γ .

The spatial cell density c(t, x) can be obtained as the marginal distribution of the structured cell density c(t, x) = Υ×P×Γ ĉ(t, x, ξ, y, α)dξ dy dα.

(

The dynamics of the structured cell density is described by

∂ ∂t ĉ(t, x, ξ, y, α) = Ŝ(t, x, ξ, y, α) -∇ x • F (t, x, ξ, y, α) -∇ ξ • Ĝ(t, x, ξ, y, α) -∇ y • Ĥ(t, x, ξ, y, α)-∇ α • K(t, x, ξ, y, α) (3)
whose full derivation is based upon work by Domschke et al. [START_REF] Domschke | Structured models of cell migration incorporating molecular binding processes[END_REF] and is given in A, along with a novel derivation of a structural source term, where Ŝ is a source term and where F , Ĝ, Ĥ, K are space-structure fluxes conjugated to the variables x, ξ, y, α, respectively.

We then proceed to more clearly define each of the flux terms in (3) as follows.

Spatial Flux

The general form of the spatial flux equation is commonly obtained from Fick's law and is given by

F (t, x, ξ, y, α) = -D c ∇ x ĉ(t, x, ξ, y, α) +ĉ(t, x, ξ, y, α)χ v ∇ x v(t, x) +ĉ(t, x, ξ, y, α) q i=1 χ i (y)∇ x m i (t, x), (4)
where the first term represents the spatial undirected diffusion of cells, the second term and third terms correspond to directed haptotactic and chemotactic cell migration, respectively.

Structural Fluxes

We consider here the dynamics of a cell population in structure space. Each cell of the population is characterised by its structure state vector s = (ξ, y, α) and by its location x ∈ D. We consider that cells in the same location follow a dynamics defined by the vector field Ψ on s ∈ Υ ×P ×Γ , with c(t, x), m(t, x), v(t, x) as parameters defining the local environment

ds dt = Ψ (s; c(t, x), m(t, x), v(t, x)). (5)
Different cells have different initial conditions at t = t 0 , whose distribution is given by ĉ(t 0 , x, s). Let s(t) = Φ t,t0 (s 0 ) be the unique solution of (5) starting from s 0 at t 0 .

Let us consider the cell sub-population located in bounded spatial V ⊂ D and structural U ⊂ Υ ×P ×Γ boxes. A population in which each cell follows (5) fulfils the continuity equation, namely

V Φt,t 0 (U ) ĉ(t, x, s) ds dx = V U ĉ(t 0 , x, s) ds dx - t t0 ∂V Φ t ,t 0 (U ) F (t , x, s) • n(x) ds dσ(x) dt + t t0 V Φ t ,t 0 (U ) Ŝ(t , x, s) ds dx dt , (6)
where Φ t,t0 (U ) is the image of U by Φ t,t0 , ∂V is the boundary of V , n(x) is the normal vector and dσ(x) is the surface measure on this boundary. Performing a change of variables in the left hand side of (6) we get

Φt,t 0 (U ) ĉ(t, x, s) ds = U ĉ(t, x, Φ t,t0 (s))J t,t0 ds , (7)
where

J t,t0 = | det dΦt,t 0 ds | is the Jacobian determinant.
Using Stokes theorem and the first fundamental theorem of calculus in (6) and further using [START_REF] Gavutis | Determination of the Two-Dimensional Interaction Rate Constants of a Cytokine Receptor Complex[END_REF] it follows

V U d dt [ĉ(t, x, Φ t,t0 (s))J t,t0 ] ds dx = - V Φt,t 0 (U ) ∇ x • F (t, x, s) ds dx + V Φt,t 0 (U ) Ŝ(t, x, s) ds dx. (8)
After changing the structure variables in the two integrals in the right hand side of ( 8) we get

d dt [ĉ(t, x, Φ t,t0 (s))J t,t0 ] = -∇ x • F (t, x, s)J t,t0 + Ŝ(t, x, s)J t,t0 . (9)
Using

1 J dJ dt = ∇ s • Ψ (s, c(t, x), m(t, x), v(t, x))
, from ( 9) we obtain the Liouville equation

∂ĉ(t, x, s) ∂t = -∇ s •(ĉ(t, x, s)Ψ (s, c(t, x), m(t, x), v(t, x)))-∇ x • F (t, x, s)+ Ŝ(t, x, s).
(10) Comparing this result to Eq.3 it follows that the structural fluxes Ĝ, Ĥ, Î are advection fluxes

Ĝ = ĉΨ ξ (ξ, y, α; c(t, x), m(t, x), v(t, x)), (11)
Ĥ = ĉΨ y (ξ, y, α; c(t, x), m(t, x), v(t, x))), (12)
Î = ĉΨ α (ξ, y, α; c(t, x), m(t, x), v(t, x)), (13)
where Ψ ξ , Ψ y , Ψ α are the components of the vector Ψ on the directions ξ, y, α, respectively.

Dynamics in receptoro-binding space

Notice that each ligand binds to the available cognate receptors. Thus, the binding rate depends on the free receptor amount ξ i -y i and is proportional to the ligand concentration

m i b i (ξ, y, m) = β i ϑ(ξ i -y i )m i , (14)
where ϑ is a function allowing to cope with the situation when binding is thresholded in the concentration of free receptors. The unbinding rate is simply proportional to the fraction of the carrying capacity of bound receptors

u i (y) = η i y i . (15)
Bound receptors are internalised with a rate

ι i (y) = k i y i . (16)
A subset of these internalised receptors are recycled. The timescale ζ -1 i of this process results from complex interactions between receptors and scaffolds inside the endosome [START_REF] Grant | Pathways and mechanisms of endocytic recycling[END_REF] and depends nonlinearly on y. Therefore, the recycling rate reads

r i (y) = ζ i (y)y i , 0 ≤ ζ i (y) ≤ k i . (17)
Receptors are synthesised by the cell with a rate p i (α, ξ) that depends on the metabolic variables α and also on actual concentration of receptors ξ and are lost by various mechanisms with a rate proportional to ξ

d i (ξ) = d i ξ i . (18)
In summary, the receptoro-binding variables of a single cell follow the differential equations

dξ dt = Ψ ξ (ξ, y, α) = P(α, ξ) -Dξ + (R(y) -I)y (19)
dy dt = Ψ y (ξ, y) = βB(ξ -y)m -(U + I)y, (20)
where P, D, R, I,β, B, U are diagonal matrices with diagonal entries

p i , d i , ζ i , k i , β i , ϑ(ξ i -y i ), η i , respectively.
It follows that the advection fluxes in receptor and binding spaces are

Ĝ(t, x, ξ, y, α) = ĉ(t, x, ξ, y, α)[P(α, ξ) -Dξ + (R(y) -I)y] (21) Ĥ(t, x, ξ, y, α) = ĉ(t, x, ξ, y, α)[βB(ξ -y)m -(U + I)y]. (22)

Dynamics in metabolic space

The part of internalisation flux that is not recycled and that escapes lysosome degradation triggers signaling and induces changes of the metabolic variables α. We use a flux-based description of these variables that considers that there are γ irreversible metabolic fluxes each one producing a different molecule. The reversible case can be simply obtained by doubling the number of variables for each reversible flux. To each one of these fluxes we associate a scalar variable 0 ≤ α i ≤ 1, meaning no production activity and maximum production activity for α i = 0 and α i = 1, respectively. In order to represent competition between fluxes we impose the condition

r i=1 α i ≤ 1. Thus α ∈ Γ , where Γ = {(α 1 , . . . , α r ) | 0 ≤ α i ≤ 1, r i=1 α i ≤ 1} is a simplex.
This description is equivalent to the space of admissible fluxes in stoichiometric and flux balance analysis of metabolic networks where α i , 1 ≤ i ≤ r represent activities of extreme pathways or currents [START_REF] Clarke | Stoichiometric network analysis[END_REF][START_REF] Schilling | Theory for the systemic definition of metabolic pathways and their use in interpreting metabolic function from a pathwayoriented perspective[END_REF]. The dynamics in the metabolic space is described phenomenologically imposing the invariance of the simplex Γ as fundamental property. A possible such choice is

dα i dt = Ψ αi (y, α) = f i (y)(1 -α i ) -µ i α i , (23)
where

µ i ≥ µ 0 > 0, f i ≤ f 0 , f 0 > 0, rf 0 < (f 0 + 1)
. The corresponding advection flux in the metabolic structure space is Î = (I 1 , . . . , I γ ) with the components

I i = ĉ[f i (y)(1 -α i ) -µ i α i ] . (24)

Spatial dynamics of diffusible ligands

Begin by denoting m := [m 1 , . . . , m p , m p+1 , . . . , m q ] T , with m j := m j (t, x), as the total vector of molecular species, where there exist q molecular species of which the first p ≤ q species are binding ligands.

Then, the spatial dynamics of all molecular species are defined by a diffusive process, and with a species specific diffusion coefficient D mj for m j (t, x). The binding ligands, within the molecular species, are removed from the population of free molecules through binding. All molecules are produced by the cellular population, in a metabolic-activity-dependent manner, and are either contributed to or detracted from by a situation specific sink or source function Θ(t, x). Therefore, denoting the q-dimensional vectors of parameters • := [• 1 , . . . , • q ] T , we obtain the relations for molecular species as

∂ m ∂t =∇ x • diag(D m)∇ x m -ε Γ P Υ (diag( βϑ(ξ -y)) m -diag(η)ȳ)ĉ dξ dy dα + Γ P Υ φα (α)(1 -m)ĉ dξ dy dα + Θ(t, x) , (25)
with φα (•) : Υ → R q defining a vector of production values for each molecular species given the cellular metabolic activity level, α; β = [β 1 , . . . , β p , 0, . . . , 0] T ; and η = [η 1 , . . . , η p , 0, . . . , 0] T ; ε is a constant converting surface to volume binding/unbinding rates.

Particularised INF-Based Model

It is necessary to first have a discussion about the context into which we shall place this model, with respect to the generalised SST framework for SARs. First of all, and for simplicity, we neglect the receptoral space and source terms in the INF case. This is due to the fact that we do not consider the creation of INF SARs, but rather their behaviour and spatial recruitment, and the change in binding in the INF case appears to be related to affinity rather than flux of the binding proteins themselves.

It should be clear, that a main concern in modeling the INF system is the numerical simplification involved in the reduction of the number of necessary dimensions under consideration. This has succeeding consequences in terms of our ability to intuit the results of the system and better understand both the SST framework, and the internal processes for communicative SARs. For this reason, we also neglect, initially, the spatial dynamics of the SAR cells and concentrate first on cell-cell communication mediated by the diffusable ligand. Now, we contextually define the binding variable, y ∈ P with v = p = 1, 0 ≤ y ≤ξ = 1, such that increasing values of y correspond to the increasing concentration of bound INF-INFAR1-INFAR2 complexes for some given (t, x, α) ∈ I ×D×Γ .

The metabolic variable, α ∈ Γ , is somewhat more complicated in biological terms since we wish to encapsulate a state of the cell under which a certain reaction is more likely to take place. In the particular case of INF, for example, we wish to understand the metabolic variable as describing a state of the cell wherein ISGs, such as IRF-7, are more frequently transcribed. Therefore, begin by describing α = 0 as a state in which ISGs are not transcribed and α = 1 as some state where ISGs are transcribed at their physiologically maximal rates. Then we understand α, itself, as encapsulating the propensity for the cell to proactively transcribe ISGs through the activity of the Jak-Stat pathway.

Within this paradigm, then, these two variables will interact in the following way. Begin by considering a scenario in which one cluster of INF SARs are stimulated by a single initial dosage of INF. The cell will bind these INF molecules and increase in binding state of the cell, y, will form the INF-INFAR1-INFAR2 complex and initiate the reactions of the Jak-Stat pathway. This will subsequently increase the cells metabolic state, α, of the cell and cause the increased production of INF. The increase in transcription of ISGs will then cause a decrease in the affinity of INFAR1 and INFAR2 for INF. This feed-back effect reported in reference will decrease the maximal value of binding state, y, that can be achieved by the unstimulated cell. This, in turn, will decrease the physiological concentration of the Jak-Stat reactions and reduce the metabolic state, α, of the cell.

Unthresholded Binding Model

The binding rates have chosen illustrative values consistently with previous models [START_REF] Domschke | Structured models of cell migration incorporating molecular binding processes[END_REF][START_REF] Trucu | Multiscale computational modelling and analysis of cancer invasion[END_REF] but with the difference that we consider here the negative feedback loop of the interferon system between the metabolic state of the cell and the binding of molecular species to the surface. In this respect, we consider binding to be non-dimensionalised and that feedback causes the maximal binding rate to decrease linearly with the metabolic state of the cell such that the range of values of y for which positive binding exists is given by y < 1 -α. Thus, we consider that the binding dynamics of molecular species to the surface, b : I ×D×P ×Γ → R, can be given by b(m, y, α)

:= β(1 -y -α)m
where β is the binding rate constant for interferon. Fig. 2: Diagram describing the simplified interferon cell-regulatory system. Unbound interferon (top left) will attempt to bind unbound receptors on the surface of the cell (circle), in accordance with its affinity for these receptors, with binding rate β. Likewise, these bound interferon-receptor complexes unbind with some rate η. The other way in which the proportion of surface bound molecules may decrease in through the internalisation of interferon receptor complexes with rate δ. The internalisation of interferon, through a complicated biological pathway, leads to a metabolic switching of the interferon-producing cell infrastructure from the default state of dormant (OFF) to active (ON), in which state the cell produces greater levels of interferon with rate φ (centre). The cell infrastructure attempts to return to the default (OFF) state with a constant rate µ 0 . In the active (ON) state the conformation of receptors, in the presence of interferon, is reduced which can be modelled through the reduction of the ability of interferon molecules to bind their receptors (i.e. β ↓, bottom right).

When we say 'binding' in this context, we actually make a generalisation of the concept of 'meaningful binding' which is to say that binding is sufficient to allow for recruitment of the secondary complex (IFNAR2) and subsequent co-phosphorylation of their protein tails.

The rate of removal of bound molecular species from the surface of the cell has a first component corresponding to unbinding and a second component corresponding to internalisation and degradation of bound receptors. Therefore, we consider that the removal of species from the cell surface, d : P → R, can be given by d(y) := (η + δ)y ,

where η gives the unbinding rate of molecules from the surface of c 1 (t, x, y, α) and δ gives the rate of cellular degradation of bound interferon.

Further, we make the assumption that the gene responsible for regulating the production of interferon has a default transcriptional state of 'off', such that the gene is not transcribed unless appropriately upregulated. Therefore we arrive at a relation for the advective rate for change in metabolic profile, µ : P × Γ → R, of the cells which is given by µ(y, α) := δy(1 -α) -µ 0 α where δy is the internalisation-degradation rate (as above) and µ 0 is the intrinsic metabolic restoration rate, the purpose of which is to restore the default metabolic position of the cell α = 0. The term (1 -α) is chosen such that the metabolic state of the cell might never exceeds a maximum value normalized to one.

Production of m with respect to the metabolic state of the cell is given by the production rate function φ : Γ → R and is assumed to be of the form

φ(α) := φ α α(α -θ α ) ,
where φ α is the rate constant for metabolic production of m and θ α is some thresholding value above which the cell become metabolically active with respect to the production of interferon, m.

For reasons that will become clear in the following subsection we call this model the unthresholded binding model, which is then written

                 ∂c1 ∂t = -∇ y • [β(1 -y -α)m -(η + δ)y] c 1 -∇ α • [δ(1 -α)y -µ 0 α] c 1 ∂c2 ∂t = 0 ∂m ∂t = ∇ x D m ∇ x m - Γ P (β(1 -y -α)m -ηy)εc 1 (t, x, y, α) dy dα P Γ φ α α(α -θ α )(1 -m)c 1 (t, x, y, α) dα dy -λmc 2 .
(26)

Thresholded Binding Model

There are several alternative interpretations of potency of a ligand for meaningful binding and signaling triggering [START_REF] Kersh | High-and Low-Potency Ligands with Similar Affinities for the TCR : The Importance of Kinetics in TCR Signaling[END_REF]. One interpretation associate this potency to the product between concentration and affinity of the ligand, suggests that ligands are detected irrespective to their quality as long as their concentration is above a threshold. Thresholds in the number of triggered receptors have been observed for immune T cells [START_REF] Viola | T cell activation determined by t cell receptor number and tunable thresholds[END_REF]. The second interpretation is based on kinetic proof-reading and suggests that a minimal binding time is needed for a given ligand to trigger signaling [START_REF] François | The Case for Absolute Ligand Discrimination : Modeling Information Processing and Decision by Immune T Cells[END_REF]. The correlation between binding time characteristics and immune cell activation is confirmed by several studies [START_REF] Kersh | High-and Low-Potency Ligands with Similar Affinities for the TCR : The Importance of Kinetics in TCR Signaling[END_REF][START_REF] Gascoigne | T-cell receptor binding kinetics in t-cell development and activation[END_REF]. Furthermore, recent dynamical studies demonstrated the phosphorylation of STAT2 to follow the formation of the complex (which is more or less instantaneous, < 1 second) by approximately 8 seconds [START_REF] Löchte | Live cell micropatterning reveals the dynamics of signaling complexes at the plasma membrane[END_REF] for complete activation. All these studies suggest the intrinsic assumption that meaningful binding requires that receptor-ligand complex to be bound for at least a minimal time τ min . In general, depending on the comparison between the timescales of meaningful complex formation and dissolution and those of activation of the signaling processes it is possible that both concentration and temporal thresholds apply to the ligand recognition. We do not aim to resolve this issue here. Because our model does not account for binding time heterogeneity, we simply replace the temporal threshold by a concentration one, considering that there is a function τ b (m) relating the concentration of ligands to the binding time. Then, for some concentration m(t, x) = θ m we have that τ b (θ m ) = τ min such that θ m gives the concentration of m sufficient for effective binding of the INFAR2 protein and INFAR1-INFAR2 complex. In order to cope with this threshold effect, we rewrite the binding flux term as b(y, α, m)

:= β(1 -y -α)(m -θ m ).
Substituting this new relation back into our model, we obtain the thresholded binding model

                 ∂c1 ∂t = -∇ y •[β(1 -y -α)(m -θ m ) -(η + δ)y] c 1 -∇ α •[δ(1 -α)y -µ 0 α] c 1 ∂c2 ∂t = 0 ∂m ∂t = ∇ x D m ∇ x m - Γ P (β(1 -y -α)m -ηy)εc 1 (t, x, y, α) dy dα P Γ φ α α(α -θ α )(1 -m)c 1 (t, x, y, α) dα dy -λmc 2 , (27)

Spatially Dynamic, Thresholded Binding Model

Finally, we consider a spatially dynamic system wherein cells are able to move through the spatial domain. We choose to endow this system with 2 primary functions of migration: (i) diffusion, by virtue of immune cells' natural inclination to motility, and (ii) chemotaxis, by virtue of immune cells' ability to actively respond to an immune-response signal as a recruitment signal.

In stating this, we therefore assume that the immune cell will interpret the presence of IFN as a response to, for example, a viral threat to the body and respond to this signal by migrating towards its origin. We further assume that even in the absence of an IFN gradient, cell Brownian motion will generate spatial fluxes leaving regions of highest cell concentration.

We thusly rewrite the system as

                       ∂c1 ∂t = ∇ x • D c1 ∇ x c 1 -∇ y • [β(1 -y -α)(m -θ m ) -(η + δ)y] c 1 -χ m ∇ x • c 1 ∇ x m -∇ α • [μ + (1 -α)y -µ -] c 1 ∂c2 ∂t = 0 ∂m ∂t = ∇ x D m ∇ x m - Υ P (β(1 -y -α) -(η + ν r δ)y)εc 1 (t, x, y, α) dy dα P Υ φ α α(α -θ α )(1 -m)c 1 (t, x, y, α) dα dy -λmc 2 , (28)

Results from Numerical Simulations

Spatially static, single-cluster results were generated by simulating ( 26), whilst multi-cluster results generate by simulating [START_REF] Samuel | Antiviral Actions of Interferons.pdf[END_REF]. Spatially-dynamics results were generated by simulating [START_REF] Schilling | Theory for the systemic definition of metabolic pathways and their use in interpreting metabolic function from a pathwayoriented perspective[END_REF]. A full description of numerical techniques used, methods, and parameters for simulating this system of equations is given in B, where parameters were used as appropriate for the simulated model.

In the following we will refer to two types of numerical simulations that differ by the type of initial condition. Single cluster simulations start with a localized cell distribution having a single maximum. Multiple cluster simulations start with an initial cell distribution having several maxima periodically positioned in space.

Spatially-static, single-cluster simulations

Single-cluster results (Fig. 3 &4) demonstrate an initial rise in average binding position, c y , of the cellular population with a concurrent rise in average metabolic position, c α . In c we also observe the rise in metabolo-binding state with a focus developing at approximately (y, α) ≈ (0.45, 0.55), with a negatively graduated non-linear ridge, and a tail between the focus and (y, α) = (0, 0).

Beyond t = 20, the average distribution in the binding space remains largely static, whilst the population continues to redistribute itself into a teardrop geometry, around the average position. This indicates, firstly, that the cell is capable of sustaining its own binding state, through production, upon initial stimulation with interferon. The formation of this geometry could be as a result of the maximal concentration of producer cells being central, and thusly producing greater levels of interferon which can be bound by the population, itself.

The distribution in the metabolic space exhibits oscillation, around its average position, for all time points t ≥ 15 (Fig. 4 c α ). This oscillation is both transverse and longitudinal, and is likely to occur as a result of the SAR-cycling between the metabolic and binding states of these cells. This demonstrates the importance of the establishment of heterogeneity within the cellular population as it acts to regulate the interferon output of the system, whilst concurrently maximising metabolic expedition from the available and bound interferon supplies. Interferon producer cells do not act in unison and, indeed, use heterogeneity to co-regulate cells within such a cluster.

The final observation that one wishes to make in the results for the singlecluster case is the visible SAR-cycle displayed within the metabolo-binding space (Fig. 3 &4 c). Regions of the solution for the cellular population appear to increase their binding state of interferon; before concurrently increasing their metabolic state and slightly decreasing their binding state; subsequently decreasing their metabolic and binding states, together; and beginning this cycle, once more. Whilst the majority of the population maintains its position within the bulk of this distribution, there exist cells (or subpopulations of the cellular population) that are affected by this feedback cycle.

Spatially-static, multi-cluster simulations

In the multi-cluster results (Fig. 5678) we observe a significant difference in the behaviour of the metabolic and binding spaces, in comparison to those of the single cluster. One observes the appearance of stable regions within the metabolic space, at high values for α; a phenomenon that we term 'metabolic trapping'. In the low affinity case, where the focal point for metabolo-binding dynamics would be lower in value, this effect is likely due to the feeding back of interferon proteins between clusters that lead the internal feedback mechanism to be ineffective at downregulating the metabolic state of the cell. In the high affinity case, this is likely to be due to the high binding and retention rates, in comparison to the unbinding rate, which causes the internalisation rate to remain high.

The binding state (Fig. 5-8 c y ), on the other hand, demonstrate oscillatory dynamics which were before characteristic of the metabolic state. Upon the establishment of stable metabolic dynamics, at high values for α, one expects that the conflict between the high rates of binding (caused by high rates of production and subsequent values for free chemical concentrations) and the feedback mechanism of the metabolic gene circuitry would cause such a behaviour. Cells will attempt to bind the high levels of interferon whilst the feedback mechanism continually acts to diminish the affinity of producer cells for interferon. One should also notice that in the low affinity case (Fig. 5 &6), as opposed to the high affinity case (Fig. 7 & 8), one observes that the signal is conveyed to the neighbouring cells. This can only be achieved through the implementation of a threshold in the binding dynamics for c(t, x, y, α) and this same threshold mediates the distance at which the signal can be conveyed.

Moreover, a simple comparative between the high affinity multi-cluster (Fig. 7 &8), low affinity multi-cluster (Fig. 5 &6), and single-cluster (Fig. 3 &4) results will show that the concentrations of interferon produced by the low affinity multi-cluster system were far in excess of those in the other two cases. This is likely as a result of the cumulative production but also as a result of the production of the two, or more, clusters feeding back the interferon to one another, causing a metabolic trapping effect. This metabolic trapping is manifest as an emergence of the population at the upper boundary of the metabolic space and retention of this position. This effect is opposed to that of the metabolo-binding SAR-cycling that one observes in the single 27) for high affinity (λ = 0.5) are given for c(t, x, y, α) in the spatio-metabolic domain (1 st row, c α ), with x on the horizontal plane and α on the vertical axis; in the spatiobinding domain (2 nd row, c y ), with x on the horizontal plane and y on the vertical axis; and for m(t, x) in space (3 rd row ), for t ∈ {0, 15, 30} respectively.

cluster case and is as a direct result of inter-cluster heterogeneity, where the promotion of the primed state in one cluster will facilitate the priming of the second, and so on.

Spatially-dynamic, multi-cluster simulations

Consider, now, the numerically generated results for the system (28), with parameters given as in above sections (B). We give the simulated solutions for the high-affinity, multi-cluster IFN case (Fig. 9 & 10), only, as the spatio-metabolobinding dynamics are similar at both high and low affinities. One immediately observes the dissolution of the discrepancy between the two species in terms of their communicative capability. The high affinity SARs are able to communicate with one another under a spatially-dynamic, chemotactic regime. In order to best understand these dynamics, one must observe them in the passage of time. The chosen initial conditions impose a stimulus on the central cluster of cells, whilst peripheral clusters are in a state of metabolic relaxation (Fig. 9). The spatial dynamics of the central cluster, at early time-points, will be mainly balanced between diffusive processes and chemotactic autoaggregation. In the peripheral clusters, however, the absence of IFN means that the spatial dynamics are mainly dictated by diffusive processes.

This diffusion in the cellular population allows some small subpopulation of cells to migrate sufficiently towards the central cluster so as to overcome the thresholding in the metabolic dynamics. Coupling this subpopulation with high affinity molecules, one achieves a fast dynamics in the binding and metabolic spaces on the perimeter of the peripheral clusters (Fig. 9, t = 20). Once these peripheral subpopulations have been potentiated to the point where they are capable of producing high affinity IFN, the cluster attains an intra-cluster supply and is capable of maintaining its own levels of IFN (Fig. 9, m(20, x)), resulting in initially peaked levels of IFN concentration at peripheral sites.

In the chemotactic simulations, one can more clearly see the elements of inter-cluster oscillation as an illustration of similar intra-cluster events. One observes an initially raised production dynamics in the central clusters (Fig. 9, t = 10); followed by fast metabolic dynamics within, and a concurrent raising of the local concentrations around, the peripheral clusters (Fig. 9, t = 20); a subsequent response from the central cluster as the peripheral clusters feedback IFN to elevate binding rates (Fig. 10, t = 30); and the resolution of this oscillatory behaviour in the establishment of a quasi-equilibrium (Fig. 10, t ≥ 40), where intra-cluster dynamics prevail but result in little macroscopic change. The initial inter-cluster heterogeneity is a necessary precursive state for the establishment of this uniformity in behavioural dynamics. Moreover, the establishment of this synchronicity between the clusters leads to another effect stemming from the chemotactic dynamic. Not only are cells capable of communicating in the chemotactic paradigm but they also self-attenuate their diffusion and auto-aggregate upon the establishment of intra-cluster activation. This may have profound implications for immunity: If, as one might intuitively predict, cells who are inclined to utilise chemotactic dynamics were attracted to the first cluster, and activated at some gradualistic pace, then the infection of the organism by a pathogen would result in the accumulation of IFN excreting cells. The decay of the spatial diversity in the cells would then lead the body to become more vulnerable to infection at novel sites, as there would no longer be IFN SARs present. If, however and as predicted by our model, we have a slight diffusive process which allows the signal to be passed but followed by auto-aggregation, then the cells would remain mostly in situ and would propogate the signal without compromising their position in the event of a further wave of infection.

Discussion

The model and framework that we have herein developed is also sufficiently general so as to be useful in cases that extend beyond the IFN system and even beyond the more general category of SARs. Generality is achieved through the biologically global forms of the binding and unbinding functions as well as the particularly general form chosen for the metabolic flow function, which describes a whole metabolic pathway in a reasoned but condensed single ODE form. Moreover, this approach to modelling the proliferation of a cell brings this model much closer to the realistic scenario, by considering the effects on the cell over entirety of the cell-cycle, rather than assuming temporal continuity in behaviour.

The single-cluster model is the first of its kind to demonstrate a qualitative biological SAR-cycling between binding and metabolic dynamics of a SAR (Fig. 3 &4). More basic, or simplistic, models may be capable of producing quantitatively similar results but could not capture the mechanistic heterogeneity within biological systems which cause them to function as they do.

Alone, this illustrates the potential for SST systems to differentially mimic biological systems to a far greater accuracy than can current modelling techniques.

In terms of the biology, this model makes two important realisations: That low affinity molecules may be necessary, for the functioning of the system, in order that the concentrations of such molecules, at long range, are sufficiently high so as to activate distant clusters of producer cells. Also, the biological system actually has two important functions of heterogeneity internally, in order to self regulate clusters and maintain sensible levels of interferon, and externally between clusters, so as to convey the activation signal of one cluster by firstly priming an initially excited cluster at a distance.

The internal heterogeneity established by clusters informs one that the ability for a cluster of SARs (specifically those for IFN) to maintain optimal levels of metabolism and reciprocal output, it is necessary for some subpopulation of cells to sacrificially reduce their levels of binding. This appears to be as a consequence of the feedback between metabolism, α, and binding, y, such that as one subpopulation rapidly increases its metabolism it will feel and subsequent inhibition of its ability to bind and will sacrifice itself such that another subpopulations may rapidly increase its binding and metabolism, due to the increased availability of local IFN. This is an important effect of intrapopular heterogeneity which we term 'subpopular metabolic quiescence', and may explain several of the inter-cellular, intra-popular oscillatory events in biology.

The latter of these two realisations recognises the importance of heterogeneity to the biological system. We demonstrate that in order that a primary cluster be primed, upon excitation, it must be allowed to be internally heterogeneous such that more active cells serve to activate less active cells whilst down regulating their own activity. This is essential for maintenance of activity levels and eventually for switching the system off. We further show that this ability for one cluster to self-activate and autoregulate is essential to maintain the long range signal and activate further clusters, at a distance. This nuancing is not possible within the simple spatial model [START_REF] Barkai | Protease helps yeast find mating partners[END_REF].

Spatially dynamic results demonstrated a breakdown in the different abilities of high and low affinity IFN to affect inter-cluster cell communication. This demonstrates that communication can be achieved either by means of reducing the barriers to the travelling molecule (affinity to consumer cells) or by cellular migration, reducing the distance between SARs themselves. In biology these dynamics may occur in environments which have more freedom for the cells to migrate and may not be achievable in many instances. In cases where migration is not possible it may be advantageous to increase production of lower affinity IFN, where high affinity IFN may be advantageous otherwise, due to the resultant increase in dynamic rate.

The biological significance of these processes are underscored by the intricate intra-and inter-cluster spatial and metabolo-binding dynamics. The major features are an intra-cluster oscillatory dynamic and a post-communicative, intra-cluster auto-aggregation which may be immunologically advantageous (depending on the paradigm considered). In the paradigm where cells are capable of migration, however, one will immediately notice that any given signal is much harder to contain or confine to a local spatial domain, which may have further implications. Therefore, local biological considerations may effect the evolutionary choice of method for communication chosen. One phenomenon, observed within the multi-dimensional model, which cannot be recreated within more simple mathematical models is that of the augmentation of metabolic state, and therefore, production in the presence of inter-cluster cooperativity. In the simple models, one has a mechanism of feedback wherein a cluster will create interferon in the presence of interferon, amplifying a given local signal. This return, however, always achieves a maximal concentration and the rate is dependent only on local interferon concentration. In the SST context, one observes that the inter-cluster supply of IFN protein between clusters actually increases the metabolic state of all involved clusters causing the productions rates to increase, concurrently. This is a qualitative result which makes a qualitative difference to the final resting state of interferon concentration.

We recognise, also, that the conveyance of the signal in the low affinity cases (Fig. 5 &6) is dependent on some thresholding parameter in the binding space, and can be justified through the biological realisation that sufficiently low quantities of chemical are insufficient to bind the receptor for long enough durations so as to cause co-phosphorylation of the internal proteins. This is a further major difference between this and previous modelling techniques, since previous modelling techniques make no comment on this phenomenon. A demonstrable advantage of this modelling framework is the ability to flag up novel biological problems, not necessarily perceptible to other state-variable frameworks.

Finally, this framework is far more approachable for the biological community, in terms of understanding. The internal and inter-cluster heterogeneity described by the SST framework is relatable to biologists in a way that is conducive to dialogue. In line with this a further explanation proffered to the thresholding problem, however, could be that there two such interferon molecules involved in this process; one of high and one of low affinity. The high affinity molecules may serve to perpetuate the activation of the considered cell, or cluster, whilst the low affinity molecule may serve to convey this signal to other producer cells. This is a theme that the authors intend to explore in a further publication.

A Derivation of the Spatio-Structuro-Temporal Model with Receptor Space

Let D ⊂ R d with d ∈ {1, 2, 3} be a bounded spatial domain, I = [0, T ] ⊂ R, with T > 0 be an arbitrary time interval. Further, let Υ ⊂ R υ with υ ∈ N characterise the available binding sites, and corresponding binding space, for given receptor ξ i for i ∈ {1, ..., υ}, which may differ in structure dependent on the molecules capable of binding each ξ i and let P ⊂ R p characterise the binding space for the biological complexes bound to these available sites. Finally, let Γ ⊂ R γ with γ ∈ N characterise the metabolic subspace of each Υ whose boundary is given by the corresponding extreme currents for the effected metabolic gene network. Herein, the space defined by Υ × P shall be referred to as the elementary-state (estate) space and will give a characterisation of the total structure of an individual element's state. Now, let the variables x ∈ D ⊂ R d represent space; ξ ∈ Υ ⊂ R υ represent receptoral state; y ∈ P ⊂ R p represent the binding state of these receptors; and α ∈ Γ ⊂ R γ represent the metabolic state of these cells. Therefore, we have also that (ξ, y) ∈ Υ × P gives the receptoro-binding state of the population, at any given point in the spatial domain, D.

Further, let U , V , W be rectangles in D, Υ ×P, and Γ respectively (i.e. U ×V ×W ⊆ D×Υ ×P ×Γ ). Then defining ĉ(t) = W V U c(t, x, (ξ, y), α) dx d(ξ, y) dα [START_REF] Schreiber | The molecular basis for functional plasticity in type I interferon signaling[END_REF] the change in c := ĉ(t, x, (ξ, y), α) per unit time in the spatio-metabolo-receptoro-binding region U ×V ×W is given by

dc(t) dt = W V U Ŝ(t, x, (ξ, y), α) dx d(ξ, y) dα - W V ∂U F (t, x, (ξ, y), α) • n(x) dσ d-1 (x) d(ξ, y) dα - W U ∂V [ Ĝ(t, x, (ξ, y), α), Ĥ(t, x, (ξ, y), α)] T • n(ξ, y) dσ υ+p-1 (ξ, y) dx dα - V U ∂W K(t, x, (ξ, y), α) • n(α) dσ γ-1 (y) dx d(ξ, y) dα (30)
where σ d-1 , σ 2r-1 , and σ γ-1 are surface measures on ∂D, ∂P, and ∂Γ , respectively. Supposing, now, that F , G, H, and J, are in the class of continuously differentiable vector fields, C 1 , one can use Stokes' Theorem to write

dc(t) dt = W V U Ŝ(t, x, (ξ, y), α) dx d(ξ, y) dα - W V U ∇x • F (t, x, (ξ, y), α) dx d(ξ, y) dα - W U V ∇ (ξ,y) • [ Ĝ(t, x, (ξ, y), α), Ĥ(t, x, (ξ, y), α)] T d(ξ, y) dx dα - U V W ∇α • K(t,
x, (ξ, y), α) dα d(ξ, y) dx [START_REF] Trucu | Multiscale computational modelling and analysis of cancer invasion[END_REF] and using Lebesgue's Dominated Convergence Theorem, one can move the time derivative within the integral for ĉ

W V U ∂ ĉ ∂t dx d(ξ, y) dα = W V U Ŝ(t, x, (ξ, y), α) dx d(ξ, y) dα - W V U ∇x • F (t, x, (ξ, y), α) dx d(ξ, y) dα - V U W ∇ (ξ,y) • Ĝ(t, x, (ξ, y), α) Ĥ(t, x, (ξ, y), α) d(ξ, y) dx dα - U V W ∇α • K(t, x, (ξ, y), α) dα d(ξ, y) dx , (32)
which can be written

R d+υ+p+γ [ ∂ ĉ ∂t ]1 U×V×W (x, (ξ, y), α) dx d(ξ, y) dα = R d+υ+p+γ [ Ŝ(t, x, (ξ, y), α)]1 U×V×W (x,
where the functions on the right-hand side describe fluxes in the cellular population density.

A.1 Derivation of a Structural Source Term

The source term accounts for cell multiplication by division. It is clear that for the source term, therefore, one must consider the full, continuous transition from mother-cell to 2 daughter-cell and, finally, back to 2 second generation mother-cell. In order to do this, we consider a normalised micro-temporal scale, τ , for which a complete cell cycle occurs on τ ∈ [0, 1). In order to formulate the source term for the system of equations we write

Likewise, these formulae are used for the calculation of the solution for the interferon molecular species, m(t, x). We further apply the population-based constraint c τ +1 1 := ĉτ+1

in order to constrain growth in the population due to the advective term under condition c(t, x, ξ, y, α) ≥ 0. We can write this in the particular case give since S(t, x, ξ, y, α) = 0 and therefore we have that there is no overall change in population. Otherwise, however, this can be achieved by stepwise accumulation and conformity.

In order to compute accurate solutions to the multi-cluster distribution arrays, we denote one individual cluster as c 1,i (t, x, y) for any i ∈ {1, . . . , k}, where k is the total number of clusters or initial distributions. Then we have that the entire cellular population distribution is defined as c 1 (t, x, y, α) := k i=1 c 1,i (t, x, y, α).

Observe that from the fundamental theorem of calculus we, therefore, have

These constraints should either leave the population c(t, x, y, α) unaltered or correct for any small instabilities arising from the long-term cumulation of O(δ 2 ) spatial advective errors, which are not adequately dealt with by the predictor-corrector methodology. We also introduce the notations

B.2 Parameters

Here we give the table of parameters for the complete, SST system:

Fig. 1 :

 1 Fig. 1: missing caption

Fig. 3 :

 3 Fig.3: Single-cluster results from simulation of model[START_REF] Pedersen | Intra-and inter-islet synchronization of metabolically driven insulin secretion[END_REF] for low affinity (λ = 0.1) are given for c(t, x, y, α) in the spatio-metabolic domain (1 st row, c α ), with x on the horizontal plane and α on the vertical axis; in the spatiobinding domain (2 nd row, c y ), with x on the horizontal plane and y on the vertical axis; in the metabolo-binding domain (3 rd row, c), with α on the horizontal axis and y on the vertical axis; and for m(t, x) in space (4 th row ), for t ∈ {0, 5, 10, 15} respectively. take off figure titles; increase fontsize for axes labels, ticks labels, legend. explain legend in the caption text.

Fig. 4 :

 4 Fig. 4: Single-cluster results from simulation of model (26) for low affinity (λ = 0.1) are given for c(t, x, y, α) in the spatio-metabolic domain (1 st row, c α ), with x on the horizontal plane and α on the vertical axis; in the spatiobinding domain (2 nd row, c y ), with x on the horizontal plane and y on the vertical axis; in the metabolo-binding domain (3 rd row, c), with α on the horizontal axis and y on the vertical axis; and for m(t, x) in space (4 th row ), for t ∈ {20, 25, 30, 35} respectively.

Fig. 5 :

 5 Fig. 5: Multi-cluster results from simulation of model (27) for low affinity (λ = 0.01) are given for c(t, x, y, α) in the spatio-metabolic domain (1 st row, c α ), with x on the horizontal plane and α on the vertical axis; in the spatiobinding domain (2 nd row, c y ), with x on the horizontal plane and y on the vertical axis; and for m(t, x) in space (3 rd row ), for t ∈ {0, 15, 30} respectively.

Fig. 6 :

 6 Fig. 6: Multi-cluster results from simulation of model (27) for low affinity (λ = 0.01) are given for c(t, x, y, α) in the spatio-metabolic domain (1 st row, c α ), with x on the horizontal plane and α on the vertical axis; in the spatiobinding domain (2 nd row, c y ), with x on the horizontal plane and y on the vertical axis; and for m(t, x) in space (3 rd row ), for t ∈ {45, 60, 75} respectively.

Fig. 7 :

 7 Fig. 7: Multi-cluster results from simulation of model (27) for high affinity (λ = 0.5) are given for c(t, x, y, α) in the spatio-metabolic domain (1 st row, c α ), with x on the horizontal plane and α on the vertical axis; in the spatiobinding domain (2 nd row, c y ), with x on the horizontal plane and y on the vertical axis; and for m(t, x) in space (3 rd row ), for t ∈ {0, 15, 30} respectively.

Fig. 8 :

 8 Fig. 8: Multi-cluster results from simulation of model (27) for high affinity (λ = 0.5) are given for c(t, x, y, α) in the spatio-metabolic domain (1 st row, c α ), with x on the horizontal plane and α on the vertical axis; in the spatiobinding domain (2 nd row, c y ), with x on the horizontal plane and y on the vertical axis; and for m(t, x) in space (3 rd row ), for t ∈ {45, 60, 75} respectively.

Fig. 9 :

 9 Fig.9: Multi-cluster results from simulation of model (28) for high affinity (λ = 0.5) are given for c(t, x, y, α) in the spatio-metabolic domain (1 st row, c α ), with x on the horizontal plane and α on the vertical axis; in the spatiobinding domain (2 nd row, c y ), with x on the horizontal plane and y on the vertical axis; and for m(t, x) in space (3 rd row ), for t ∈ {0, 10, 20} respectively.

Fig. 10 :

 10 Fig. 10: Multi-cluster results from simulation of model (28) for high affinity (λ = 0.5) are given for c(t, x, y, α) in the spatio-metabolic domain (1 st row, c α ), with x on the horizontal plane and α on the vertical axis; in the spatio-binding domain (2 nd row, c y ), with x on the horizontal plane and y on the vertical axis; and for m(t, x) in space (3 rd row ), for t ∈ {30, 40, 50} respectively.

  (ξ, y), α) dx d(ξ, y) dα -R d+υ+p+γ

[

  ∇x • F (t, x, (ξ, y), α)]1 U×V×W (x, (ξ, y), α) dx d(ξ, y) dα -R d+υ+p+γ ∇ (ξ,y) • Ĝ(t, x, (ξ, y), α) Ĥ(t, x, (ξ, y), α) 1 U×V×W (x, (ξ, y), α) dx d(ξ, y) α -R d+υ+p+γ [∇α • K(t,x, (ξ, y), α)]1 U×V×W (x, (ξ, y), α) dx d(ξ, y) dα .

( 33 )

 33 There, since we have that {U ×V ×W | U, V, W -compact with piecewise smooth boundaries} is a family of generators for the Borelian σ-algebra on U ×V ×W we can denote 1 A as the indicator function for any arbitrary A ⊆ D×Υ ×P ×Γ and writeR d+υ+p+γ [ ∂ ĉ ∂t ]1 A (x, (ξ, y), α) dx d(ξ, y) dα = R d+υ+p+γ [ Ŝ(t, x, (ξ, y), α)]1 A (x, (ξ, y), α) dx d(ξ, y) dα -R d+υ+p+γ [∇x • F (t, x, (ξ, y), α)]1 A (x, (ξ, y), α) dx d(ξ, y) dα -R d+υ+p+γ ∇ (ξ,y) • Ĝ(t, x, (ξ, y), α) Ĥ(t, x, (ξ, y), α) 1 A (x, (ξ, y), α) dx d(ξ, y) dα -R d+υ+p+γ [∇α • K(t, x, (ξ, y), α)]1 A (x,(ξ, y), α) dx d(ξ, y) dα (34) for any arbitrary Borelian set A in the σ-algebra on D×Υ ×P ×Γ . Then we can replace 1 A with any simple function, as so R d+2r+γ [ ∂ ĉ ∂t ]ν(x, (ξ, y), α) dx d(ξ, y) dα = R d+υ+p+γ [ Ŝ(t, x, (ξ, y), α)]ν(x, (ξ, y), α) dx d(ξ, y) dα -R d+υ+p+γ [∇x • F (t, x, (ξ, y), α)]ν(x, (ξ, y), α) dx d(ξ, y) dα -R d+υ+p+γ ∇ (ξ,y) • Ĝ(t, x, (ξ, y), α) Ĥ(t, x, (ξ, y), α) ν(x, (ξ, y), α) dx d(ξ, y) dα -R d+υ+p+γ

[

  ∇α • K(t, x, (ξ, y), α)]ν(x, (ξ, y), α) dx d(ξ, y) dα∀ν ∈ C ∞ 0 (D×P ×Γ ) .(35)Then, since this relation holds for any C ∞ test function, ν(x, (ξ, y), α), we obtain the equation∂ ĉ ∂t = Ŝ(t, x, (ξ, y), α) -∇x • F (t, x, (ξ, y), α) -∇ (ξ,y) • [ Ĝ(t,x, (ξ, y), α), Ĥ(t, x, (ξ, y), α)] T -∇α • K(t, x, (ξ, y), α) ,

WŜ 2 [0, 1 )( 2 - 2 [0, 1 )( 2 -

 212212 (t, x, y) dy = 2 [0,1) (2-τ )W φ(ỹ, c, v)ĉ(t, x, ỹ) dỹ dτ -W φ(y, c, v)ĉ(t, x, y) dy (37) and use the change of variable ỹ(y) = (2 -τ )y dỹ = (2 -τ ) dy (38) such that we can write W Ŝ(t, x, y) dy = τ ) p W φ((2 -τ )y, c, v)ĉ(t, x, (2 -τ )y) dy dτ -W φ(y, c, v)ĉ(t, x, y) dy = W τ ) p φ((2 -τ )y, c, v)ĉ(t, x, (2 -τ )y) dτ dy -W φ(y, c, v)ĉ(t, x, y) dy .

∇x • c 1

 1 (t, x, y, α)∇xm = ∇x • k i=1 c 1,i (t, x, y, α) ∇xm = k i=1 ∇x • c 1,i (t, x, y, α)∇xm (45)Since we have that the overall population does not change with respect to changes time (S(t, x, y, α) = 0), we can use that (44) and (45) imply that the population constraint holds on each individual cluster of interferon producer cells

  , x, y, α) dy and cy := Γ c(t, x, y, α) dα as quantifying the spatio-metabolic and spatio-binding distributions, respectively, and c := D c(t, x, y, α) dx as quantifying the non-spatial metabolo-binding distribution of the cellular population c(t, x, y, α).

Table 1 :

 1 Table of parameters

Acknowledgements References

A.2 Derivation of a Structural Source Term for Systems with Receptors

In order to formulate the source term for the system of equations we write

and use the change of variable

such that we can write

Differentiating both sides of this relation, we arrive at our final expression of source flux for the total population as

B Numerical Methods & Parameters

B.1 Numerical methods

We use the 4 th order Runge-Kutta predictor for this system, given by cτ+1

where F (c τ 1 ) := F (c τ 1 , m τ ) are given by the local central difference approximation of the spatio-structural dynamics for c τ 1 := c 1 (t τ , x, y, α) at the given time point t τ . We then use a MacCormack corrector, of the form ) .