Figure 4.32 (page 169):

Species cccDNA versus time for hepatitis B virus model; deterministic and average stochastic models.

Code for Figure 4.32

Text of the GNU GPL.

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(:,5)]);
axis ([0,200,0,25]);
title ('Figure 4.32')

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);

Matlab dynamically-linked MEX functions (C language)

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;
}

Octave dynamically-linked functions (C++ language)

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;
}