All the reproduced figures are plotted using COPASI.
Fig 5A with 5A.sedml and model.xml xi_A=0.006 and beta_3=0.35x10^(-9)
Fig 5B with 5B.sedml and model.xml xi_A=0.003 and beta_3=0.35x10^(-9)
Fig 5C with 5C.sedml and model.xml xi_A=0 and beta_3=0.35x10^(-9)
Fig 5D with 5D.sedml and model.xml xi_A=0 and beta_3=0.28x10^(-9)

Professor Fassoni has confirmed us with some parameters used to reproduce the figures:
Because the different orders of magnitude involved( large populations and small parameters), to simulate it, we used the normalized (nondimensional) populations. 

So in the figure we plotted, we used :n(t)=N(t)/K, G(t)=G(t)/K, a(t)=A(t)/K and n(t),g(t),a(t) follow the same ODEs of N(t),G(t),A(t) (system 2.4). 

We have the following replacements in parameters: 
- r_N is replaced by r_N/K
- K_A is replaced by K_A/K
- \xi is replaced by \xi/K
- \beta_i is replaced by \beta_i * K (i=1,2,3)
- the other parameters remain the same
- the initial conditions, which are N(0)=(r_N/\mu_N)-1, G(0)=1, A(0)=0, are replaced by n(0)=N(0)/K, g(0)=G(0)/K, a(0)=A(0)/K

We used the value K=10^8.
Correction for parameter K_A in table 1 of the paper: K_A=10^7