diff --git a/N803_model_2.m b/N803_model_2.m
index 1632fee..a176570 100644
--- a/N803_model_2.m
+++ b/N803_model_2.m
@@ -1,16 +1,14 @@
 %% N803_model_2.m - solves ODE model for N-803 treatment of SIV
 %
 %     /--------------------------------------------------------------\
-%     | Date:         11/05/2021                                     |
+%     | Date:         11/27/2021                                     |
 %     | Author:       Jonathan Cody                                  |
 %     | Affiliation:  Purdue University                              |
 %     |               Weldon School of Biomedical Engineering        |
 %     |               Pienaar Computational Systems Pharmacology Lab |
 %     \--------------------------------------------------------------/
 %
-% Version 2D+ (also allowing regulation calibration)
-%
-% 12/03/2021 - added some extra terms to the full output matrix
+% Version 2D++ (replacing mu with fA)
 %
 % Nomenclature: V  = SIV virions [#/µL]
 %               E0 = resting SIV-specific CD8+ T cells [#/µL]
@@ -64,7 +62,7 @@
 %	FUNCTION
 %  ========================================================================
 function Y_OUT = N803_model_2(SoluTimes,DoseTimes,Pars,...
-                              EquatedPars,fullOut,timeOut,warnOff)
+                              EquatedPars,fullOut,timeOut,warnOff,validate)
 
 % convert inputs to column vectors
 if size(SoluTimes ,1)==1 ; SoluTimes = SoluTimes' ; end
@@ -91,7 +89,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)
-mu   = Pars(04) ;% initial E+B turnover relative to healthy NHP
+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]
@@ -141,8 +139,19 @@
 % 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),[-8 2]) ;
+x = fzero(@(x)icFun(x),[-20 10]) ;
 [~,aEi,aBi,Ei,Bi] = icFun(x) ;
 
 % regulation initial value
@@ -254,18 +263,6 @@
 % 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
 %  ========================================================================
@@ -285,8 +282,8 @@
     Ei = Sol(1:9) ;
     Bi = Sol(10:12) ;
     
-    % ensure initial conditions satisfy mu
-    targ = d*(Ei(1)+Bi(1)) + dA*(Ei(9)+Bi(3)) - mu*d*EBi ;
+    % ensure initial conditions satisfy fA
+    targ = Ei(1) + Bi(1) - (1-fA)*EBi ;
 end
 
 function dY = system(~,Y)