Figure 9.25:

Envelope versus time for hepatitis B virus model; initial guess and estimated parameters fit to data.

Code for Figure 9.25

Text of the GNU GPL.

main.m


  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
% jbr, elh,  11/28/01
%
% A = cccDNA, B = rcDNA, C = env
%
% Reactions
% 1: A -> A + B
% 2: B -> A
% 3: A -> A + C
% 4: A -> 0 (degradation)
% 5: C -> 0 (degradation)
% 6: B + C -> 0 (secreted virus)

% converted to use parest, jbr, 4/24/2007
% converted to use paresto (casadi), jbr, 4/13/2018
%

more off

model=struct;
model.nlp_solver_options.ipopt.linear_solver = 'ma27';
%model.nlp_solver_options.ipopt.print_level = 0;
%model.nlp_solver_options.print_time = false;
%model.transcription = 'shooting';

model.x = {'ca', 'cb', 'cc'};
model.p = {'k1', 'k2', 'k3', 'k4', 'k5', 'k6', 'wa', 'wb', 'wc'};
model.d = {'ma', 'mb', 'mc'};



model.ode = @hbv_rxs;
model.lsq = @(t, y, p) {p.wa*(y.ca-y.ma), p.wb*(y.cb-y.mb), p.wc*(y.cc-y.mc)};

% Set the reaction rate constant vector kr; use log transformation
kr_ac = [2; 0.025; 1000; 0.25; 1.9985; 7.5E-6];
logkr = log10(kr_ac);
p.k1 = logkr(1); p.k2=logkr(2); p.k3 = logkr(3);
p.k4 = logkr(4); p.k5 = logkr(5); p.k6=logkr(6);

% output times
tfinal = 100;
ndata = 51;
tdata = linspace(0, tfinal, ndata)';
model.tout = tdata;

% measurement weights
mweight = sqrt([1; 1e-2; 1e-4]);
p.wa = mweight(1); p.wb = mweight(2); p.wc = mweight(3);
p_ac = [logkr; mweight];

pe = paresto(model);

% Set the initial condition
small = 0;
x0_ac  = [1; small; small];

y_ac = pe.simulate(zeros(3, 1), x0_ac, p_ac);

rng(1);
R = diag([0.1 0.1 0.1])^2;
noise = sqrt(R)*randn(size(y_ac));
% proportional error
y_noisy = y_ac .* (1 + noise);
y_noisy = max(y_noisy, 0);

% list of all parameters
theta0 = [logkr; mweight; x0_ac];
lbtheta = theta0;
ubtheta = theta0;

% estimate kr(1)-kr(6); loosen bounds
est_ind = 1:6;
loose = 5.5;
lbtheta(est_ind) = -loose*ones(numel(est_ind),1);
ubtheta(est_ind) = loose*ones(numel(est_ind),1);

more off
[est, y, p] = pe.optimize(y_noisy, theta0, lbtheta, ubtheta);

theta_conf = pe.confidence(est, est_ind, 0.95);

disp('Optimal parameters and confidence intervals')
[est.theta(est_ind), theta_conf]

% % optimal fit
% figure(1)
% plot(tdata, est.x, tdata, y_noisy, 'o')

% figure(2)
% semilogy(tdata, est.x, tdata, y_noisy, 'o')

% % diagnose problems

% r_ac = diag(mweight)*(y_ac-y_noisy);
% phi_ac = sum(diag(r_ac'*r_ac))
% rfit = diag(mweight)*(est.x-y_noisy);
% phi_fit = sum(diag(rfit'*rfit))

% check eigenvalues/eigenvectors for badly determined parameter
% directions

[u, lam] = eig(est.d2f_dtheta2(est_ind,est_ind))

%
%  watch out, eigenvalues may be unordered;
%  run sort so i(1) will be the index of the smallest eigenvalue;
%  still have a +/- sign ambiguity on u(:,i(1)), but thetap with either
%  sign on u(:,i(1)) should give the same fit to the data

[s,i] = sort(diag(lam));
logkrp = est.theta(est_ind) + 0.5*u(:,i(1));

% try text, page 545
%logkrp = [0.32; -1.47; 4.21; -0.46; 1.56; -5.08];

thetap = theta0;
thetap(est_ind) = logkrp;
lbtheta = thetap;
ubtheta = thetap;

more off
[estp, yp, pp] = pe.optimize(y_noisy, thetap, lbtheta, ubtheta);
% simulate with parameters perturbed in weak direction

store1 = [model.tout, [y.ca;y.cb;y.cc]', est.x', estp.x'];
store2 = [model.tout, y_noisy'];
save hbv_det.dat store2 store1;

if (~ strcmp (getenv ('OMIT_PLOTS'), 'true')) % PLOTTING
subplot (3, 2, 1);
plot (store2(:,1), store2(:,2), '+', store1(:,1), store1(:,[2,5]));
axis ([0, 100, 0, 60]);
% TITLE hbv_det_cccdna

subplot (3, 2, 3);
plot (store2(:,1), store2(:,3), '+', store1(:,1), store1(:,[3,6]));
axis ([0, 100, 0, 600]);
% TITLE hbv_det_rcdna

subplot (3, 2, 5);
plot (store2(:,1), store2(:,4), '+', store1(:,1), store1(:,[4,7]));
axis ([0, 100, 0, 30000]);
% TITLE hbv_det_env

subplot (3, 2, 2);
plot (store1(:,1), store1(:,[2,8]));
axis ([0, 100, 0, 60]);
% TITLE hbv_det_cccdnap

subplot (3, 2, 4);
plot (store1(:,1), store1(:,[3,9]));
axis ([0, 100, 0, 600]);
% TITLE hbv_det_rcndap

subplot (3, 2, 6);
plot (store1(:,1), store1(:,[4,10]));
axis ([0, 100, 0, 30000]);
% TITLE hbv_det_envp
end % PLOTTING

hbv_rxs.m


1
2
3
4
5
function rhs = hbv_rxs(t, y, p)
  kr = 10.^([p.k1, p.k2, p.k3, p.k4, p.k5, p.k6]);
  rhs = {kr(2)*y.cb - kr(4)*y.ca, ...
         kr(1)*y.ca - kr(2)*y.cb - kr(6)*y.cb*y.cc, ...
         kr(3)*y.ca - kr(5)*y.cc - kr(6)*y.cb*y.cc};