%% mydata_Venturia_canaescens: parasitic wasp
% Ana L Llandres (04/05/2013)

%% References
% * <http://www.bio.vu.nl/thb/deb/deblab/add_my_pet/add_my_pet.pdf *add_my_pet manual*> 
% * <http://www.bio.vu.nl/thb/research/bib/LikaKear2011.html *LikaKear2011*>
%   Lika, K., Kearney, M. R., Freitas, V., Veer, H. W. van der, Meer, J. van der, Wijsman, J. W. M., Pecquerie, L. and Kooijman, S. A. L. M.
%   The `covariation method' for estimating the parameters of the standard Dynamic Energy Budget model {I}: philosophy and approach.   
%   J. Sea Res. 66: 270--277.
% * <http://www.bio.vu.nl/thb/research/bib/Kooy2010.html *Kooy2010*>
%   Kooijman, S.A.L.M. 2010 
%   Dynamic Energy Budget theory for metabolic organisation. Cambridge Univ. Press
%   Table 8.1, page 300
%  * Harvey, J.A., Harvey, I.F. and Thompson, D.J. 1994. Flexible Larval Growth Allows Use of a Range of Host Sizes by a Parasitoid Wasp. Ecology 75: 1420-1428.
%  * Harvey, J. A., Harvey, I.F. and Thompson, D.J. 2001. Lifetime
%    reproductive success in the solitary endoparasitoid, Venturia canescens. Journal of Insect Behavior 14: 573-593. 
%  * Ohbayashi, T., Iwabuchi, K. and Mitsuhashi, J. 1994. In Vitro rearing
%    of a larval endoparasitoid, Venturia canescens (Gravenhorst)(Hymenoptera: Ichneumonidae). I Embryonic Development. Appl Entomol Zool, 29: 123-126.
%  * Howell, J. and Fisher, R.C. 1977. Food conversion efficiency of a parasitic wasp Nemeritis canescens. Ecological Entomology 2: 143-151. 
%  * Spanoudis, C. G., and S. S. Andreadis. 2012. Temperature-dependent survival, development, and adult longevity of the koinobiont endoparasitoid 
%    Venturia canescens (Hymenoptera: Ichneumonidae) parasitizing Plodia interpunctella (Lepidoptera: Pyralidae). Journal of Pest Science 85:75-80.


%% Facts
% * Venturia canescens is a parasitoid wasp species that present four life stages: embryo (egg), juvenile (larvae), pupa and adult (wasp) stages

COMPLETE = 2.5;   % judge the level using LikaKear2011; adjust if you have less or more data
FIT = 8.9;        % compute after having obtained the estimates
                      
  global dwm % pass d_O, w_O, mu_O directly to predict_Venturia_canescens

% set data

% zero-variate data
% real data

ab = 4;          %  1 d, age at hatching of larvae when feeding at host instar 5 (age 0 is at onset of embryo development) (Harvey et al. 1994)
  T_ab = 273 + 25; % K, temperature for ah (Harvey et al. 1994)
aj = 10;         %  2 d, age at pupation when feeding at host instar 5 (Harvey et al. 1994)
  T_aj = 273 + 25; % K, temperature for aj (Harvey et al. 1994)
ae = 21.53;      %  3 d, age at emergence of adult from the pupa when feeding at host instar 5 (Harvey et al. 1994)
  T_ae = 273 + 25; % K, temperature for ap (Harvey et al. 1994)
  %am= 17;          % 4 d, given adult average lifespan (eliminating the
  %post-reproductive lifespan)
tm= 36.5;          % 4 d, given adult average lifespan of 39 emerged wasps that fed in honey solution add libitum (37.9 days, Harvey et al. 1994), this is an estimation of life span when feeding at host instar 5 (37.9d +21.53d)  T_am = 273 + 25; % K, temperature for am (Harvey et al. 1994)
    T_tm = 273 +25;
W0 = 1.87e-007;    % 5 g, dry weight of the recently oviposited egg, calculated from the lenght (l= 0.027 cm, ) and width (w= 0.0047 cm radio=r2 0.0047/2= 0.0024) of a egg picture
                    %and by assuming that the egg is cylindrical and the
                    %density of the eggs is 0.4 g/cm^3, Mass = density*pi*l*r2^2 =0.4 *pi*0.027* 0.0024^2= 1.87e-007g )
                    %* pi * 0.0047^2/4 (Ohbayashi et al. 1994)
Wb = 6.9e-5;    %  6 g, dry weight at hatching at Host Instar5  (Harvey et al. 1994)
Wj = 4.77e-3;   %  7 g, dry weight at pupation at Host Instar5 (Harvey et al. 1994)
We = 1.75e-3;   %  8 g, dry weight at emergence at Host Instar5 (Harvey et al. 1994)
del_S = 0.19;       %  9 -, mg cocoon/mg pupa =1.063 /5.495= 0.19, proportional mg (dry weight) of cocoon from the total mg of the pupa including the cocoon. Caclulated from data publisehd by Howell & Fisher 1977. (table 1).  (Howell and Fisher 1977)
N_e= 40.204;      %  10 -, number of eggs at emergence (40.2041) (Harvey et al. 2001)

% pseudo-data from pars_my_pet at T_ref from Kooy2010
v = 0.02;     % 12 cm/d, energy conductance
kap = 0.8;    % 13-, allocation fraction to soma = growth + somatic maintenance
kap_R = 0.95; % 14 -, reproduction efficiency
p_M = 18;     % 15 J/d.cm^3, [p_M] vol-specific somatic maintenance
p_T =  0;     % 16 J/d.cm^2, {p_T} surface-specific som maintenance
k_J = 0.002;  % 17 1/d, maturity maintenance rate coefficient
kap_G = .8;   % 18 -, growth efficiency

% pack data
data = [ab; aj; ae; tm; W0; Wb; Wj; We; del_S; N_e;  % 01:10 real data
        v; kap; kap_R; p_M; p_T; k_J; kap_G];                % 11:17 pseudo data

% weight coefficients for WLS and ML criteria differ
Data = [data(:,[1 1]), min(100,1 ./ max(1e-6, data) .^ 2)]; % nmregr, nrregr (WLS criterion)
Data(1:10,3) = 10 * Data(1:10,3); % give real data more weight
Data(9:3) = 200 * Data(9:3);
Data(5:8,3) = 10 * Data(5:8,3); % give weight data more weight
Data(17,3) = 200 * Data(17,3);  % more weight to kap_G

% insert temperature data for rates and times in the first column
Data([1 2 3 4], 1) = [T_ab; T_aj; T_ae; T_tm];


txt_data = {... % for presentation of predictions
    	'1 ab, d, age at hatching when fed on host instar 5 ';
        '2 aj, d, age at pupation when fed on host instar 5' ;
        '3 ae, d, age at emergence ' ;
    	'4 tm, d, life span ';
        '5 W0, g, dry weight of recently laid egg';
        '6 Wb, g, dry weight at hatching when fed on host instar 5';
    	'7 Wj, g, dry weight at pupation when fed on host instar 5';
    	'8 We, g, dry weight at puberty when fed on host instar 5';
    	'9 del_S, -, weight cocoon shell/weight pupa (included cocoon)';
        '10 N_e,-, number eggs at emergence';
        '11 v, cm/d, energy conductance ';
        '12 kap, -, allocation fraction to soma ';
        '13 kap_R, -, reproduction efficiency ';
    	'14 [p_M], J/d.cm^3, vol-spec som maint ';
    	'15 {p_T}, J/d.cm^2, sur-spec som maint ';
        '16 k_J, 1/d, maturity maint rate coefficient ';
        '17 kap_G, -, growth efficiency'};


% uni-variate data at f = 0.8 and T = 273 + 25

% larval growth (Harvey et al. 1994, graph 4 larvae growth when feeding in instar 5 of the host)
tWd_bj =  [... % time since hatching (d), larval dry weight (mg)
0	0.069444444 
1	0.069444444
2	0.222222222
3	0.722222222
4	3.180555556
5	3.666666667
6	4.777777778];
tWd_bj(:,2) = tWd_bj(:,2)/1e3; % convert mg to g
tWd_bj= tWd_bj(:,[1 2 2]);
tWd_bj(:,3) = 1000*1e-1/mean(tWd_bj(:,2))^2;
% pupal growth (Harvey et al. 1994 graph 4 pupa growth when after feeding in instar 5 of the host)
tWd_je =  [... % time since pupation (d), pupal dry weight (mg)
0	4.777777778
1	4.069444444
2	3.972222222
3	2.777777778
4	2.361111111
5	2.736111111
6	2.541666667
7	3.083333333
8	2.805555556
9	2.541666667];
tWd_je(:,2) = tWd_je(:,2)/1e3; % convert mg to g
tWd_je = tWd_je(:,[1 2 2]);
tWd_je(:,3) = 10000*1e-1/mean(tWd_je(:,2))^2;

% Harvey et al. 2001 
% data represent averages of 10 inidividuals each day
% data collected when animals face high food quantity (costant food access to 50% honey solution)and 24 hours
%   of daily host access until wasp death. Note that this species do not host-feed, so the host access is just to oviposit.
tN_em=  [... % time since adult emergence (d), reprod rate (number progeny/d), 
1	40.20408163
2	30.40816327
3	30.20408163
4	33.36734694
5	32.14285714
6	28.46938776
7	23.7755102
8	18.97959184
9	9.183673469
10	2.857142857
11	3.87755102
12	0.918367347
13	0.408163265
14	0
15	0.306122449
16	0.714285714
17	0.816326531
18	0
19	0
20	0.204081633
21	0
22	0
23	0
24	0
25	0.204081633
26	0
27	0
28	0
29	0
30	0
31	0
32	0
33	0
34	0];
tN_em(:,2)=tN_em(:,2)*100/83; % 83% of the eggs survive until adulthood which means that we need to multiply the progeny by 100/83 to get the real number of eggs(Spanoudis & Andreadis 2012)
tN_em(:,2) = cumsum(tN_em(:,2)); % cumulative reproduction
tN_em = [tN_em 0 * tN_em(:,2) + 1000/mean(tN_em(:,2))^2];


% conversion coefficients from Kooy2010 
% chemical indices  from Kooy2010 Fig 4.15
%       X     V     E     P
n_O = [1.00, 1.00, 1.00, 1.00;  % C/C, equals 1 by definition
       1.80, 1.80, 1.80, 1.80;  % H/C  these values show that we consider dry-mass
       0.50, 0.50, 0.50, 0.50;  % O/C
       0.15, 0.15, 0.15, 0.15]; % N/C
%       C     H     O     N
n_M = [ 1     0     0     0;    % C/C, equals 0 or 1
        0     2     0     3;    % H/C
        2     1     2     0;    % O/C
        0     0     0     1];   % N/C
    
% specific densities
%       X     V     E     P
d_O1 = [0.4;  0.4;  0.4;  0.4];    % g/cm^3, specific densities for organics in the embryo stage (before feeding)

d_O = [0.2;  0.2;  0.2;  0.2];    % g/cm^3, specific densities for organics

% chemical potentials from Kooy2010 Tab 4.2
%        X     V     E     P
mu_O = [525; 500;  550;  480] * 1000; % J/mol, chemical potentials for organics

% molecular weights
w_O = n_O' * [12; 1; 16; 14];  % g/mol, mol-weights for organics

% pack coefficients
dwm = [d_O, w_O, mu_O, d_O1]; % g/cm^3, g/mol, kJ/mol spec density, mol weight, chem pot



% parameters: initial values at T_ref zhen f=0.95
T_ref  = 293;      %  1 K, temp for which rate pars are given; don't change this vulue
T_A  = 8000;       %  2 K, Arrhenius temp 
f = 0.95;             %  3 -, scaled functional response
z =3.979 ;        %  4 -, zoom factor; for z = 1: L_m = 1 cm
del_M =  0.1520999026;  %  5 -, shape coefficient
F_m = 6.5;         %  6 l/d.cm^2, {F_m}, max spec searching rate
kap_X = 0.94;      %  7 -, digestion efficiency of food to reserve
v =   0.0207 ;       %  8 cm/d, energy conductance
kap =   0.9309;     %  9 -, allocation fraction to soma = growth + somatic maintenance
kap_R = 0.85;      % 10 -, reproduction efficiency
p_M =     245.7  ;       % 11 J/d.cm^3, [p_M], vol-specific somatic maintenance
p_T =  0;          % 12 J/d.cm^2, {p_T}, surface-specific som maintenance
k_J = 0.002;       % 13 1/d, maturity maint rate coefficient
E_G =  2.121e+004 ;        % 14 J/cm^3, [E_G], spec cost for structure
E_Hx=    0.0002938    ;    % 15 J,maturity at pre-birth (when embryo starts feeding)
E_Hb =   0.04483     ;    % 16 J, E_H^b, maturity at birth
y_EV_l = 0.8;      % 17 mol/mol, yield of imago reserve on larval structure
E_He   =  5.566 ;     % 18 J, E_H^p, maturity at emergence
%E_H0= 5.566;   % 19 J, maturity at start of egg maturation
E_H0 = 5.48;   % 19 J, maturity at start of egg maturation
%h_a =  0.041;      % 20 1/d^2, Weibull aging acceleration
h_a =  0.0045;      % 20 1/d^2, Weibull aging acceleration
s_G = 1e-4;        % 21 -, Gompertz stress coefficient
s_j = 0.9831;        % 22 -, reprod buffer/structure at pupation as fraction of max
s_s = 0.1057;         % 23 -, reprod buffer/structure at silk production, as fraction of max [ER]
kap_SA = 0.03;     % 24 -, weight coefficient for assimilation
kap_SG = 0.5;     % 25 -, weight coefficient for assimilation
kap_SD = 0.94;     % 26 -, weight coefficient for assimilation
k_E_im = 0.2911;   % 27 1/d, reserve turnover and reproduction buffer turnover in the imago
% pack parameters and fix T_ref and f and possibly other at well
%   in second column: 0 = fix; 1 = release

    
pars = [T_ref 0; T_A 0; f    0; z      1; del_M 0; F_m  0;  
        kap_X 0; v   0; kap  1; kap_R  0; p_M   1; p_T  0;  
        k_J   0; E_G 0; E_Hx 1; E_Hb 01; y_EV_l 0; E_He  0; E_H0 0; h_a  0; 
        s_G   0; s_j 0; s_s 0; kap_SA 0; kap_SG 0; kap_SD 0; k_E_im 0];


txt_pars = { ...    % for presentation of parameter estimates
  'T_ref, K'; 'T_A, K';          'f, -'; 
  'z, -';     'del_M, -';        '{F_m}, l/d.cm^2'; 
  'kap_X, -'; 'v, cm/d';         'kap, -'; 
  'kap_R, -'; '[p_M], J/d.cm^3'; '{p_T}, J/d.cm^2'; 
  'k_J, 1/d'; '[E_G], J/cm^3'; 'E_Hx, J';   'E_Hb, J'; 
  'y_EV_l, mol/mol';'E_He, J'; 'E_H0, j'; 'h_a, 1/d^2'; 's_G, -'; 
  's_j, -'; 's_s, -'; 'kap_SA, -'; 'kap_SG, -'; 'kap_SD, -'; 'k_E_im, 1/d'};

% estimate parameters, see LikaKear2011
nmregr_options('default'); % set options for parameter estimation
nmregr_options('max_step_number',5e3); % set options for parameter estimation
nmregr_options('max_fun_evals',2e4);   % set options for parameter estimation

%pars = nmregr('predict_Venturia_canescens', pars, Data, tWd_bj, tWd_je, tN_em);    % WLS estimate parameters using overwrite
sd = 0 * pars(:,1);                                  % initiate standard deviations

%[cov cor sd] = pregr('predict_Venturia_canescens', pars, Data, tWd_bj, tWd_je, tN_em); % get standard deviation for WLS

% get FIT
Data(:,3) = 0; Data(1:10,3) = 1; % give unit weight to real data, zero to pseudo-data
[MRE RE] = mre('predict_Venturia_canescens', pars, Data, tWd_bj, tWd_je, tN_em); % WLS-method
FIT = 10 * (1 - MRE); % get mark for goodness of fit 

% get predictions
tbj = linspace(0,6,800)'; % times for plotting lavae data
tje = linspace(0,9,800)'; % times for plotting pupa data
tem = linspace (1,34, 800)'; % times for plotting reproduction data
[Edata EW_bj EW_je EN_em] = predict_Venturia_canescens(pars(:,1), Data, tbj, tje, tem); % notice use of first column of pars only

%% Results

fprintf('FIT = %3.1f \n', FIT)
fprintf('COMPLETE = %3.1f \n\n', COMPLETE)

printpar(txt_data, data, Edata, 'data and expectations'); % for zero-variate data
fprintf('\n') % insert a blank line in screen output
printpar(txt_pars, pars, sd);

close all % remove existing figures, else you get more and more if you retry

figure % one figure to show results of uni-variate data

subplot(3,1,1)
plot(tWd_bj(:,1), tWd_bj(:,2), '.b', tbj, EW_bj, 'b')
xlabel('time since larvae hatching, d')
ylabel('dry weight, g')


subplot(3,1,2)
plot(tWd_je(:,1), tWd_je(:,2), '.b', tje, EW_je, 'b')
xlabel('time since pupation, d')
ylabel('dry weight, g')

subplot(3,1,3)
plot(tN_em(:,1), tN_em(:,2), '.b', tem, EN_em, 'b')
xlabel('time since wasp emergence, d')
ylabel('cummulated number of progeny, individuals')
