Computational Methods to Model Persistence

Abstract

Bacterial persister cells are dormant cells, tolerant to multiple antibiotics, that are involved in several chronic infections. Toxin–antitoxin modules play a significant role in the generation of such persister cells. Toxin–antitoxin modules are small genetic elements, omnipresent in the genomes of bacteria, which code for an intracellular toxin and its neutralizing antitoxin. In the past decade, mathematical modeling has become an important tool to study the regulation of toxin–antitoxin modules and their relation to the emergence of persister cells. Here, we provide an overview of several numerical methods to simulate toxin–antitoxin modules. We cover both deterministic modeling using ordinary differential equations and stochastic modeling using stochastic differential equations and the Gillespie method. Several characteristics of toxin–antitoxin modules such as protein production and degradation, negative autoregulation through DNA binding, toxin–antitoxin complex formation and conditional cooperativity are gradually integrated in these models. Finally, by including growth rate modulation, we link toxin–antitoxin module expression to the generation of persister cells.

Appendices

Appendix 1: Numerical Code to Solve an ODE/SDE

A simple matlab code to solve Eq. 10, with (D ≠ 0) or without noise (D = 0), can be found here below:

1  function ODE_SDE

2 

3 %% parameters

4 prodAT    = 0.0530 * 0.116086/0.00203;

5 degrAT    = 2.8881e-4;

6 D               = 25;

7 dt             = 0.01;        % [s] simulation time step

8 dt_save  = 10;            % [s] plotting time step

9 t_end      = 5*60*60;  % [s] final time

10 

11 %% initialize system

12 AT              = 0;

13 t_saved   = [];

14 AT_saved = [];

15 count       = 0;

16 

17 %% simulate the stochastic differential equation

18 for n = 0:((t_end)/dt)

19        t = n * dt;

20 

21     %% Euler-Heun

22  noise = sqrt(D) * sqrt(dt) * randn(); % sample from the noise

23 

24     AT_star = AT + dt * F(AT) + noise;

25     AT      = AT + (dt/2)*(F(AT) + F(AT_star)) + noise/2;

26 AT = max(AT, 0); % force protein concentration to be positive

27 

28     %% Save data

29         if (count == dt_save/dt)

30         t_saved(end+1) = t;

31         AT_saved(end+1) = AT;

32         count = 0;

33     end

34     count = count + 1;

35 

36 end

37 

38 %% plot the results

39 figure;

40 plot(t_saved./3600, AT_saved, 'k');

41 xlabel('Time (h)')

42 ylabel('AT')

43 

44 %% definition of the differential equation

45     function dATdt = F(TA)

46         dATdt = prodAT - degrAT*TA;

47     end

48 end

Appendix 2: Numerical Code Using the Gillespie Algorithm

A simple matlab code to solve Eq. 1 using the stochastic Gillespie Algorithm can be found here below:

1  % Gillespie code

2 % There are 2 reactions and there is one species AT

3 

4 %% Parameters

5 prodAT = 0.0530*0.116086/0.00203; % reaction 0  -> AT

6 degrAT = 2.8881e-4;                               % reaction AT -> 0

7 

8 %% Initialization

9 AT                  = 0;              % [AT] initial concentration AT

10 t                  = 0;              % [s] starting time

11 t_end         = 5*60*60; % [s] final time

12 t_saved    = [];            % [s] stored times

13 AT_saved   = [];

14 

15 %% Simulation

16 while t <= t_end

17     %% Update propensities

18     p1 = degrAT * AT;

19     p2 = prodAT;

20 

21     %% Computation of the random time step

22     p0  = p1 + p2;

23     r1  = rand();

24     r2  = rand();

25     dt  = 1/p0 * log(1/r1); % [s] next time step

26 

27     %% Selection of random reaction

28     %% Update the population based on selected reaction

29     yr2 = r2 * p0;

30     if yr2 <= p1

31         % reaction 1

32         AT = AT - 1;

33     else

34         % reaction 2

35         AT = AT + 1;

36     end

37 

38     %% Update the current time

39     t = t + dt;

40 

41     %% Save population information

42     t_saved(end+1) = t;

43     AT_saved(end+1) = AT;

44 

45 end

46 

47 %% plot the results

48 figure;

49 plot(t_saved./3600, AT_saved, 'k');

50 xlabel('Time (h)')

51 ylabel('AT')