From 69a766384bec7e8e2154af388240912464f2618c Mon Sep 17 00:00:00 2001
From: "Cody, Jonathan William" <cody2@purdue.edu>
Date: Mon, 24 Jan 2022 14:20:37 -0500
Subject: [PATCH] Add files via upload

---
 N803_model_2.m | 36 ++++++++++++++++++++----------------
 1 file changed, 20 insertions(+), 16 deletions(-)

diff --git a/N803_model_2.m b/N803_model_2.m
index a176570..08819e3 100644
--- a/N803_model_2.m
+++ b/N803_model_2.m
@@ -1,14 +1,17 @@
 %% N803_model_2.m - solves ODE model for N-803 treatment of SIV
 %
 %     /--------------------------------------------------------------\
-%     | Date:         11/27/2021                                     |
+%     | Date:         11/05/2021                                     |
 %     | Author:       Jonathan Cody                                  |
 %     | Affiliation:  Purdue University                              |
 %     |               Weldon School of Biomedical Engineering        |
 %     |               Pienaar Computational Systems Pharmacology Lab |
 %     \--------------------------------------------------------------/
 %
-% Version 2D++ (replacing mu with fA)
+% Version 2D+ (also allowing regulation calibration)
+%
+% 12/03/2021 - added some extra terms to the full output matrix
+% 11/27/2021 - replaced mu with fA -> initial frequency: (Ea+Ba)/(E+B)
 %
 % Nomenclature: V  = SIV virions [#/µL]
 %               E0 = resting SIV-specific CD8+ T cells [#/µL]
@@ -62,7 +65,7 @@
 %	FUNCTION
 %  ========================================================================
 function Y_OUT = N803_model_2(SoluTimes,DoseTimes,Pars,...
-                              EquatedPars,fullOut,timeOut,warnOff,validate)
+                              EquatedPars,fullOut,timeOut,warnOff)
 
 % convert inputs to column vectors
 if size(SoluTimes ,1)==1 ; SoluTimes = SoluTimes' ; end
@@ -89,7 +92,7 @@
 Vi   = Pars(01) ;% V initial value [log(#/mL)]
 EBi  = Pars(02) ;% E+B initial value [#/µL]
 fE   = Pars(03) ;% initial frequency: E/(E+B)
-fA   = Pars(04) ;% initial frequency (Ea+Ba)/(E+B)
+fA   = Pars(04) ;% initial frequency: (Ea+Ba)/(E+B)
 q    = Pars(05) ;% V growth rate (if E+B were absent) [/d]
 psi  = Pars(06) ;% Ba/Ea killing rate ratio [gB0/gE0]
 EB50 = Pars(07) ;% 50% E+B proliferation saturation [#/µL]
@@ -139,19 +142,8 @@
 % ratio between initial activation rate of E and of B
 phi = (kB*mB/(mB+dA) - 1) / (kE*mE/(mE+dA) - 1) ;
 
-% !!!!! VALIDATION --------------------------------------------------------
-if ~isempty(validate)
-    Xs = linspace(-20,20,1000) ;
-    TARG = zeros(1000,1) ;
-    for n1 = 1:1000
-        TARG(n1) = icFun(Xs(n1)) ;
-    end
-    plot(Xs,TARG,'k') ; grid on
-end
-% !!!!! VALIDATION --------------------------------------------------------
-
 % solve for aE,aB initial values and E and B initial values
-x = fzero(@(x)icFun(x),[-20 10]) ;
+x = fzero(@(x)icFun(x),[-15 5]) ;
 [~,aEi,aBi,Ei,Bi] = icFun(x) ;
 
 % regulation initial value
@@ -263,6 +255,18 @@
 % Ba killing proportion
 Y_OUT(:,37) = gB0*Ba./(1+gB2*R) ./ Kill ;
 
+% Average E proliferation rate
+Y_OUT(:,38) = ( p0*F(:,15).*E0 + pE*sum( Y(:,03:9) ,2) ) ./ (E0+Ea) ;
+
+% Average B proliferation rate
+Y_OUT(:,39) = ( p0*F(:,15).*B0 + pB*Y(:,12) ) ./ (B0+Ba) ;
+
+% Average E killing rate
+Y_OUT(:,40) =  gE0*Ea./(1+gE2*R) ./ (E0+Ea) ;
+
+% Average B killing rate
+Y_OUT(:,41) =  gB0*Ba./(1+gB2*R) ./ (B0+Ba) ;
+
 %% ========================================================================
 %	NESTED FUNCTIONS
 %  ========================================================================