## Code for Figure 4.35

### main.m

```%% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt
%%
%% This program is free software; you can redistribute it and/or
%% modify it under the terms of the GNU General Public License as
%% published by the Free Software Foundation; either version 2, or (at
%% your option) any later version.
%%
%% This program is distributed in the hope that it will be useful, but
%% WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
%% General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program; see the file COPYING.  If not, write to
%% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
%% MA 02111-1307, USA.
%%
%%
%% Hepatitis B model, stochastic simulation

%% Hepatitis B model, stochastic simulation.  Data and functions defined
%% in this file are used in both hbv_stoch.m and hbv_binary.m

global k
k = [1; 0.025; 1000; 0.25; 2; 7.5e-6];
stoi = [0, 1, 0; 1, -1, 0; 0, 0, 1; -1, 0, 0; 0, 0, -1; 0, -1, -1];
x0  = [1; 0; 0];
stoiT = stoi';
tfin  = 200;
nts   = 200;
tout = linspace(0,tfin,nts)';
%%
%% steady state analysis
%%
cas = ( k(1)-k(4) ) /   ...
( k(6) *k(4)/k(2)*( (-k(1)+k(4) +k(3)) / (k(5))) );
cbs = k(4)/k(2)* cas;
ccs = ( k(1)-k(4) ) / ( k(6) *k(4)/k(2) );
Jac = [-k(4), k(2), 0;
k(1), -k(2)-k(6)*ccs, -k(6)*cbs;
k(3), -k(6)*ccs, -k(5)-k(6)*cbs];
eig(Jac);

%%
%% deterministic dynamics
%%
opts = odeset ('AbsTol', sqrt (eps), 'RelTol', sqrt (eps));
[tsolver, x] = ode15s(@hbv_deriv, tout, x0, opts);

%  function xout = stochsim(x0, k, stoiT, tout)
%    nts  = length(tout);
%    tfin = tout(nts);
%    time = tout(1);
%    x = x0;
%    iout = 1;
%    while (time < tfin)
%      r(1) = k(1)*x(1);
%      r(2) = k(2)*x(2);
%      r(3) = k(3)*x(1);
%      r(4) = k(4)*x(1);
%      r(5) = k(5)*x(3);
%      r(6) = k(6)*x(2)*x(3);
%      rtot = sum(r);
%      p=rand(2,1);
%      %% choose likely time of next reaction
%      tau = -log(p(1))/rtot;
%      timenew = time + tau;
%      %% choose which reaction (mth) is likely to occur
%      rcum = 0;
%      m = 0;
%      while ( rcum <= p(2)*rtot)
%        m = m + 1;
%        rcum = rcum + r(m);
%      end
%      xnew  = x + stoiT(:,m);
%      %% check if passing time of requested output
%      while( tout(iout) <= timenew )
%        xout(:,iout) = x;
%        iout = iout + 1;
%        if (iout > nts) break; end
%      end
%      time = timenew;
%      x = xnew;
%    end
%  endfunction
%xplot1 = stochsim(x0, k, stoiT, tplot);
%%
%% check the deterministic steady state and eigenvalues
%%
%% check parameters for hbv_det simulations
%%k(1) = 2;
%%k(5) = 1.9985

%%
%% already ran simulation, just load file and write a plotting file
%%

load hbv_binary.mat;

%%
%% save for plotting
%%
%% nullseed(116) generates a nice figure
%%
rand('seed', nullseed(116))
xnull = stochsim(x0, k, stoiT, tout);

table = [tout, x, xavg', xout', xnull'];

plot (table(:,1),[table(:,2),table(:,8),table(:,11)]);
axis ([0,200,-4,35]);
title ('Figure 4.35')
```

### hbv_deriv.m

```function rhs = hbv_deriv(t,x);
global k
rhs = zeros (3, 1);
rhs(1) = k(2)*x(2) - k(4)*x(1);
rhs(2) = k(1)*x(1) - k(2)*x(2) - k(6)*x(2)*x(3);
rhs(3) = k(3)*x(1) - k(5)*x(3) - k(6)*x(2)*x(3);
```

### stochsim.c

```/*

Copyright (C) 2001, James B. Rawlings

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License as
published by the Free Software Foundation; either version 2, or (at
your option) any later version.

This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; see the file COPYING.  If not, write to
the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA.

*/

#include
#include

#include

static double
get_rand (void)
{
double rnd;

mxArray *plhs[1];

mexCallMATLAB (1, plhs, 0, 0, "rand");

rnd = (mxGetPr (plhs[0]))[0];

mxDestroyArray (plhs[0]);

return rnd;
}

void
mexFunction (int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[])
{
int nargin = nrhs;

const mxArray *xin;    double *pxin;
const mxArray *k;      double *pk;
const mxArray *stoiT;  double *pstoiT;  int stoiT_nr;
const mxArray *tout;   double *ptout;

mxArray *xout;   double *pxout;

int nts, nx;
double tfin, time;
int iout = 0;
double x0, x1, x2;
double k0, k1, k2, k3, k4, k5;
int pxout_offset = 0;

if (nargin < 4 || nargin > 4)
mexErrMsgTxt ("usage: stochsim (x0, k, stoiT, tout)");

/* Maybe we should check that these are numeric, to avoid crashing
Matlab if someone passes something else?  */
xin = prhs[0];
k = prhs[1];
stoiT = prhs[2];
tout = prhs[3];

pxin = mxGetPr (xin);
pk = mxGetPr (k);
pstoiT = mxGetPr (stoiT);
ptout = mxGetPr (tout);

stoiT_nr = mxGetM (stoiT);

nts = mxGetNumberOfElements (tout);
tfin = ptout[nts-1];

time = ptout[0];

nx = mxGetNumberOfElements (xin);

if (nx != 3)
mexErrMsgTxt ("stochsim: expecting X to be a 3-element vector");

x0 = pxin[0];
x1 = pxin[1];
x2 = pxin[2];

xout = mxCreateDoubleMatrix (nx, nts, mxREAL);
pxout = mxGetPr (xout);

k0 = pk[0];
k1 = pk[1];
k2 = pk[2];
k3 = pk[3];
k4 = pk[4];
k5 = pk[5];

while (time < tfin)
{
double r0 = k0 * x0;
double r1 = k1 * x1;
double r2 = k2 * x0;
double r3 = k3 * x0;
double r4 = k4 * x2;
double r5 = k5 * x1 * x2;

double rtot = r0 + r1 + r2 + r3 + r4 + r5;

double timenew, rnd, rcum, tmp;
int m = 0;
double xnew0, xnew1, xnew2;

/* return if system has extinguished */
if (rtot == 0.0)
{
while (iout < nts)
{
pxout[pxout_offset++] = x0;
pxout[pxout_offset++] = x1;
pxout[pxout_offset++] = x2;

iout++;
}
break;
}

/* choose likely time of next reaction */
timenew = time - log (get_rand ()) / rtot;

/* choose which reaction (mth) is likely to occur */
rnd = get_rand ();
rcum = r0;
m = 0;
tmp = rnd*rtot;
if (rcum <= tmp)
{
++m; rcum += r1;

if (rcum <= tmp)
{
++m; rcum += r2;

if (rcum <= tmp)
{
++m; rcum += r3;

if (rcum <= tmp)
{
++m; rcum += r4;

if (rcum <= tmp)
++m; rcum += r5;
}
}
}
}

xnew0 = x0 + pstoiT[stoiT_nr*m];
xnew1 = x1 + pstoiT[stoiT_nr*m+1];
xnew2 = x2 + pstoiT[stoiT_nr*m+2];

/* check if passing time of requested output */
while (ptout[iout] <= timenew)
{
pxout[pxout_offset++] = x0;
pxout[pxout_offset++] = x1;
pxout[pxout_offset++] = x2;

iout++;

if (iout >= nts)
break;
}

time = timenew;

x0 = xnew0;
x1 = xnew1;
x2 = xnew2;
}

plhs[0] = xout;
}
```

### stochsim.cc

```/*

Copyright (C) 2001, James B. Rawlings

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License as
published by the Free Software Foundation; either version 2, or (at
your option) any later version.

This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; see the file COPYING.  If not, write to
the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA.

*/

#include
#include
#include
#include
#include
#include

DEFUN_DLD (stochsim, args, ,
"stochsim(x0, k, stoiT, tout)")

{
octave_value retval;

int nargin  (args.length ());

if (nargin < 4 || nargin > 4)
{
print_usage ("stochsim");
return retval;
}

const ColumnVector xin (args(0) . vector_value ());
const ColumnVector k (args(1) . vector_value ());
const Matrix stoiT (args(2) . matrix_value ());
const ColumnVector tout (args(3) . vector_value ());

if (error_state)
{
error ("stochsim: invalid arguments");
return retval;
}

const int nts  = tout.length ();
const double tfin = tout(nts-1);

double time = tout(0);
int iout = 0;

const int nx = xin.length ();

if (nx != 3)
{
error ("stochsim: expecting X to be a 3-element vector");
return retval;
}

double x0 = xin(0);
double x1 = xin(1);
double x2 = xin(2);

Matrix xout (nx, nts, 0.0);

double *pxout = xout.fortran_vec ();

int pxout_offset = 0;

double k0 = k(0);
double k1 = k(1);
double k2 = k(2);
double k3 = k(3);
double k4 = k(4);
double k5 = k(5);

while (time < tfin)
{
// Allow this loop to be interrupted.
OCTAVE_QUIT;

double r0 = k0 * x0;
double r1 = k1 * x1;
double r2 = k2 * x0;
double r3 = k3 * x0;
double r4 = k4 * x2;
double r5 = k5 * x1 * x2;

double rtot = r0 + r1 + r2 + r3 + r4 + r5;

// return if system has extinguished
if (rtot == 0.0)
{
while (iout < nts)
{
pxout[pxout_offset++] = x0;
pxout[pxout_offset++] = x1;
pxout[pxout_offset++] = x2;

iout++;
}
break;
}

// choose likely time of next reaction
double timenew = time - log (octave_rand::scalar ()) / rtot;

// choose which reaction (mth) is likely to occur
double rnd = octave_rand::scalar ();
double rcum = r0;
int m = 0;
double tmp = rnd*rtot;
if (rcum <= tmp)
{
++m; rcum += r1;

if (rcum <= tmp)
{
++m; rcum += r2;

if (rcum <= tmp)
{
++m; rcum += r3;

if (rcum <= tmp)
{
++m; rcum += r4;

if (rcum <= tmp)
++m; rcum += r5;
}
}
}
}

double xnew0 = x0 + stoiT(0,m);
double xnew1 = x1 + stoiT(1,m);
double xnew2 = x2 + stoiT(2,m);

// check if passing time of requested output
while (tout(iout) <= timenew)
{
pxout[pxout_offset++] = x0;
pxout[pxout_offset++] = x1;
pxout[pxout_offset++] = x2;

iout++;

if (iout >= nts)
break;
}

time = timenew;

x0 = xnew0;
x1 = xnew1;
x2 = xnew2;
}

retval = xout;

return retval;
}
```