Mathematics and Computational Methods for Complex Systems summative-assessment.pdf Prof. Luc Berthouze

Introduction

This coursework is worth 100% and is due Thursday, 6 January 2022, 16:00 (but

please always check what Canvas and Sussex Direct says!). Submission is electronic, via Canvas. Submission guidelines are identical to the formative coursework,

namely: (a) Analytical work can be hand-written or typeset; (b) simulation code +

analysis can be presented separately or integrated into a Jupyter / Colab notebook

(whatever you prefer!). Unlike the previous coursework, there is no page limit but

conciseness of exposition is always best!

In this coursework, you are going to engage in some proper modelling work. You are

given a scenario along with a couple of ordinary differential equations that describe

the operation of this system (in principle I could have asked you to derive them from

the scenario description but I don’t). The questions take you through a number of

steps by which you will gain a detailed understanding of the system / model. Although

all questions are related, it is possible to get results even if you find yourself unable

to complete some of the steps.

Statement of problem

We consider a system of individuals (the nature of which is irrelevant) that can take

two states A and B. At any point in time, an individual can only be in one state.

There are no other states possible. We will denote by [A] the number of individuals

in state A and by [B] the number of individuals in state B. There are N individuals

in total. All individuals can be assumed to be in contact with all other individuals. At

time t = 0, there are B0 individuals in state B. All other individuals are therefore in

state A. Individuals in state A turn into state B at a rate β [B]
N , i.e., proportional to the

number of their neighbours that are in state B. Individuals in state B turn into state A

at a fixed rate γ.

Now, in an ideal world, you should be able to write the two ordinary differential equations (ODE) that describe this system. To make sure that everyone can proceed with

the coursework, I am providing the ODEs in the last page of this document. Before

consulting that section, you are strongly encouraged to try and come up with your

own formulation.

1

Mathematics and Computational Methods for Complex Systems summative-assessment.pdf Prof. Luc Berthouze

Questions

We start with some analytical work.

Analytical work

1. Using the fact that [A] + [B] = N at all times, write ˙ [B] as a function of [B], i.e.,

the expression should no longer involve [A]. This is your (so-called) mean-field

equation.

(2 marks)

2. Find the equilibria of the system and determine their stability. From now on,

we will refer to the non-zero equilibrium as B∗. You may find it useful to write

your results in terms of the following quantity R0 = β

γ . Plot the phase portrait,

i.e., [B] vs [A], identifying the equilibria and their stability (following convention

described in the 2nd synchronous lecture of Unit 5).

(8 marks)

3. Produce the bifurcation plot for this system, that is, plot the value of the equilibria

as a function of R0, with R0 taking values from 0.1 to 5.0. For this question, the

value of N is irrelevant (provided it’s strictly positive) so use 1000 for example.

This should be done using Python.

(4 marks)

4. In this question, you are going to integrate ˙ [B] analytically to obtain an expression for [B](t), i.e., an expression that gives the number of individuals in state

B in time. It is rarely the case that this can be done but with this system, it is

possible. You will do this in four steps:

• Starting from the mean-field equation, factorise the right hand side by [B]
2,

then write an expression for 1

[B]2 ˙ [B]. (4 marks)

• Consider the following variable substitution: y = 1

[B]
. Using the chain rule,

express y˙ in terms of [B], then derive a simple expression for y˙, i.e., this expression should only involve y terms and parameters of the system. There

shouldn’t be any [B] or [A]. However, it will be helpful to use B∗ (calculated

in the 2nd question) to simplify the expression. (4 marks)

2

Mathematics and Computational Methods for Complex Systems summative-assessment.pdf Prof. Luc Berthouze

• Integrate this equation. You should be able to do this without any help, but

if help is needed, you should note that this expression looks very much like

the equation we solved during a synchronous lecture in Unit 4, replacing λ

and I by appropriate quantities. You can then use the result to derive an

equation for y(t). Please see short document summarising the derivation

from the lecture. (4 marks)

• You can now produce a fully worked out expression for [B](t) by remembering that [B] = B0 at time t = 0. (4 marks)

NB: In going through this question, do make sure to consider all scenarios possible regarding the value of R0.

5. Using different values of B0 (between 1 and N – briefly discuss the case B0 = 0),

plot solutions of [B](t) for various values of R0 between 0.1 and 5.0 (with γ = 0.5

for example). Confirm your expression for [B](t) is correct by (a) verifying that

it converges to B∗ for large times t and (b) visually confirming agreement when

integrating the mean-field equation using Euler (use Python). What happens

when R0 = 1? Speculate as to what this means. We will get back to this. For a

given value of R0, what happens when the value of γ changes? Provide a brief

explanation.

(11 marks)

Simulation work

Some more could be done analytically but let’s now turn to simulations. For this

part of the work, we will employ the Gillespie algorithm, an algorithm often used to

simulate complex systems. There is something very important to understand here.

The mean-field equation provided to you is derived from considering the interactions

of a very large number of individuals. The Gillespie algorithm, instead, provides

discrete simulations of the system with few individuals by explicitly simulating every

transitions. In other words, a single run of the Gillespie algorithm represents one

sample trajectory of all possible trajectories for this system. In principle, the average

of a (sufficiently) large number of Gillespie runs should converge to the mean-field.

This is what you will test here. NB: Gillespie provides a mathematically rigorous and

computationally efficient alternative to agent-based modelling. You are not asked to

compare the two approaches even if that would be an interesting thing to do (NB: but

certainly not within this assessment!).

To make sure everybody can work, I provide a Python implementation of the Gillespie

algorithm for this problem. Use the code although feel free to write your own if you feel

3

Mathematics and Computational Methods for Complex Systems summative-assessment.pdf Prof. Luc Berthouze

so inclined – you would need to read up on the Gillespie algorithm obviously. Function

gillespie_ABA takes 5 arguments: N, B0, β, γ and Tmax (max simulation time) and

produces three outputs: T, A and B which contain the time of each transition (when

one, and only one, individual changed state, whatever that state was), the numbers

of individuals in state A and in state B after that transition, respectively. NB: The

Gillespie algorithm uses exponentially-distributed times to next event so each run will

have a different timeline T of events. This won’t be apparent to you if plotting with

lines but will matter in what comes below.

1. Explore the behaviour of the system when considering suitably chosen scenarios, i.e., focus on the limit cases (e.g., small R0, large R0 and R0 = 1; small

N, large N; small B0, large B0). For each scenario, use the code provided to

generate many realisations of the stochastic process. Plot all realisations on a

single plot. Make relevant qualitative observations.

(10 marks)

2. For each scenario, calculate the average (and standard deviation) of the realisations. Here, you are going to face a problem linked with the nota bene from

the introductory paragraph. You will need to think of a solution and implement

it. Superimpose the average (and error bars) to the realisations. Use a larger

line width for visibility.

(11 marks)

3. Finally, superimpose the mean-field solution. Again, use a larger line width and

different colour for visibility. Describe and interpret agreement between average

of stochastic realisations and mean-field in relation to the choice of parameters.

In this question, using B0 = 1 (i.e., only one individual in state B at t = 0) can

help exacerbate the differences and help you think about what is happening.

You may want to refer to your bifurcation plot.

(8 marks)

4. [Slightly challenging question]: Consider 100 replications for N = 1000, β =

0.51, γ = 0.5 and 100 replications for N = 1000, β = 0.95, γ = 0.5. You should

notice a substantial difference in agreement between the mean-field and the average of the stochastic realisations depending on which scenario is considered.

How could you improve agreement for the scenario with the poorest agreement.

Please note: The difference in B∗ is not the quantity of interest here. Rather

you should think about why agreement is so poor. This does not actually involve

analytical work. An excellent answer would see you implement your proposed

solution and provide evidence of improved agreement.

(12 marks)

4

Mathematics and Computational Methods for Complex Systems summative-assessment.pdf Prof. Luc Berthouze

A bit of critical thinking

In this last part of the coursework, you are invited to think a bit more about what you

have just done.

1. So far the brief has provided no information whatsoever regarding what the

states A and B are and what the individuals are. Thinking about what is happening in this system, provide at least one example of real-world scenario to

which this model could apply. Bonus points will be given for any answer that

provides two examples, one in which the equilibria are of interest and the other

in which the critical regime (when R0 = 1) is of interest. Either way, what is the

benefit of being able to model the system?

(8 marks)

2. [Very tough question]: The model provided implicitly assumes that all individuals are potentially in contact with each other. What would be a more likely

scenario? What changes would have to be made to the code of the Gillespie

algorithm in order to include such a scenario? If you are able to do this, do

it. Then, speculate as to what could affect the results observed in the previous

questions. If you feel so inclined, demonstrate it experimentally. NB: Only 10

marks have been given to this question. However, anyone managing it successfully would receive an extra 10 marks for the assignment (with the total

mark capped to 100 obviously).

(10 marks)

Total marks: 100.

5

Mathematics and Computational Methods for Complex Systems summative-assessment.pdf Prof. Luc Berthouze

The ODEs

The system is fully described by the following two ordinary differential equations:

˙ [B] = β [B]
N [A] − γ[B]
˙ [A] = −β [B]
N [A] + γ[B]
where β and γ are strictly positive real numbers. The first equation provides the rate

of change of the number of individuals in state B and says that it is the sum of the

number of individuals in state A that become B with a rate β [B]
N (proportional to the

number of other individuals in state B) minus the number of individuals in state B that

turn to state A with rate γ. I leave it to you to interpret the second equation. Summing

both equations lead to 0 which is consistent with the fact that the total population

[A] + [B] = N at all times. To complete the description we need the initial conditions.

We will call B0 the number of individuals in state B at t = 0. The number of individuals

in state A is therefore N − B0.

6