diff --git a/N803_model_1.m b/N803_model_1.m
index 2ba0e1b..39659b5 100644
--- a/N803_model_1.m
+++ b/N803_model_1.m
@@ -2,7 +2,7 @@
 % and NK cells during a subcutaneously administered regimen of N803
 
 %     /--------------------------------------------------------------\
-%     | Date:         06/29/2019                                     |
+%     | Date:         12/05/2019                                     |
 %     | Author:       Jonathan Cody                                  |
 %     | Affiliation:  Pienaar Computational Systems Pharmacology Lab |
 %     |               Weldon School of Biomedical Engineering        |
@@ -56,7 +56,7 @@
 
 %     P(21) = killing stimulation factor []
 %     P(22)
-%     P(23) = proliferation stimulation factor []
+%     P(23) = maximum expansion rate [/d] (will convert to rho)
 %     P(24)
 %     P(25) = tolerance effect factor []
 %     P(26)
@@ -65,6 +65,7 @@
 %     P(28)
 %     P(29) = proliferation regulation factor []
 %     P(30)
+%     P(31) = tolerance recovery []
 
 %% ========================================================================
 %	OPTIONAL INPUTS
@@ -100,6 +101,8 @@
 function [Y_MAIN,extra] = N803_model_1(SoluTimes,DoseTimes,P,...
                                        EquatedPars,timeOut,normOut)
 
+warning off
+                                   
 % convert inputs to column vectors
 if size(SoluTimes ,1)==1 ; SoluTimes = SoluTimes' ; end
 if size(DoseTimes ,1)==1 ; DoseTimes = DoseTimes' ; end
@@ -109,9 +112,11 @@
 % times, and end time (used for piecewise solution of ODEs)
 PieceIn = [ 1 ; find(ismember(SoluTimes,DoseTimes)) ; length(SoluTimes) ] ;
 
-% declare index in Y (see 'system' function) for C, last tolerance
-% variable, and last regulation variable
-TolRegIn = [ 6 ; 6+P(10) ; 6+P(10)+P(11) ] ;
+% declare index in Y (see 'system' function) for tolerance variables
+TolIn = 7:(6+P(10)) ;
+
+% declare index in Y (see 'system' function) for regulation variables
+RegIn = (7+P(10)):(6+P(10)+P(11)) ;
 
 % equate parameters in 'P' based on indices in 'EquatedPars'
 if ~isempty(EquatedPars) ; P = P(EquatedPars) ; end 
@@ -126,8 +131,11 @@
 P(01) = 10^(P(01)-3) ;% [log(#/mL)] to [#/µL]
 P(16) = P(16)*P(15)  ;% [*P(15)] to [#/µL-d]
 
+% convert maximum expansion rate [/d] to prolif. stimulation factor []
+P(23:24) = P(23:24)./P(17:18) ;
+
 % declare initial conditions for Y (see 'system' function)
-Yic      = zeros(1,TolRegIn(3))                ;% zero except for...
+Yic      = zeros(1,6+P(10)+P(11))              ;% zero except for...
 Yic(1:4) = [ P(01)*[1-P(02) P(02)] P(03:04)' ] ;% V,W,E,K
 
 % declare proliferation rate constants 'q_V','q_W','r_E','r_K' [/d]
@@ -171,10 +179,8 @@
     % Y(4) = K = NK cells [#/µL]
     % Y(5) = X = N803 at absorption site [pmol/kg]
     % Y(6) = C = N803 plasma concentration [pM]
-    % Y(7      :5+P(10)      ) = tolerance delay variables
-    % Y(        6+P(10)      ) = tolerance acting variable
-    % Y(7+P(10):5+P(10)+P(11)) = regulation delay variables
-    % Y(        6+P(10)+P(11)) = regulation acting variable
+    % Y( 7 : (6+P(10)) ) = tolerance variables
+    % Y( (7+P(10)) : (6+P(10)+P(11)) ) = regulation variables
     
     % throw error if 'solvertime' exceeds 'timeOut' limit
     if toc(solvertime) >= timeOut
@@ -200,17 +206,24 @@
     %% --------------------------------------------------------------------
     % Calculate rates of change for tolerance & regulation variables ------
     
-    % for tolerance and regulation
-    for m = 1:2
-        % if number of variables is not zero
-        if P(09+m) ~= 0
-            % update variable rates of change
-            dY(TolRegIn(m)+1:TolRegIn(m+1)) = P(11+m) * ...
-                ( [ Y_terms(1) ; Y(TolRegIn(m)+1:TolRegIn(m+1)-1) ] ...
-                -                Y(TolRegIn(m)+1:TolRegIn(m+1))     ) ;
-        end
+    % if number of tolerance variables is not zero
+    if P(10) ~= 0
+        % update variable rates of change
+        dY(TolIn(1:end-2)) = P(12) * [ Y_terms(1) ; Y(TolIn(1:end-3)) ] ...
+                           - P(12) *                Y(TolIn(1:end-2)) 	;
+        dY(TolIn(end-1)) = P(12) *             Y(TolIn(end-2)) ...
+                         - P(12) * (1+P(31)) * Y(TolIn(end-1)) ;
+        dY(TolIn(end)) = P(12) * P(31) * Y(TolIn(end-1)) ...
+                       - P(12) * P(31) * Y(TolIn(end))   ;
     end
-    
+ 
+    % if number of regulation variables is not zero
+    if P(11) ~= 0
+        % update variable rates of change
+        dY(RegIn) = P(13) * [ Y_terms(1) ; Y(RegIn(1:end-1)) ] ...
+                  - P(13) *                Y(RegIn)            ;
+    end
+
 end
 
 %% ------------------------------------------------------------------------
@@ -219,7 +232,7 @@
 function Y_terms = terms(Y)
     
     % percolate vector for storing terms
-    Y_terms = zeros(12,1) ;
+    Y_terms = zeros(11,1) ;
     
     % Y_terms(1) = N803 effect []
     Y_terms(1) = Y(6) / ( P(09) + Y(6) ) ;
@@ -227,14 +240,14 @@
     % Y_terms(2) = N803 effect for E (with tolerance) []
     % Y_terms(3) = N803 effect for K (with tolerance) []
     Y_terms(2:3) = [Y_terms(1);Y_terms(1)]./ ...
-                   ( 1 + P(25:26) * Y(TolRegIn(2)) ) ;
+                   ( 1 + P(25:26) * sum(Y((5:6)+P(10))) ) ;
     
     % Y_terms(4) = killing modifier for E (with effect & regulation) []
     % Y_terms(5) = killing modifier for K (with effect & regulation) []
     % Y_terms(6) = proliferation modifier for E (w/ effect & regulation) []
     % Y_terms(7) = proliferation modifier for K (w/ effect & regulation) []
     Y_terms(4:7) = ( 1 + P(21:24).* [Y_terms(2:3);Y_terms(2:3)] )./ ...
-                   ( 1 + P(27:30) * Y(TolRegIn(3)) ) ;
+                   ( 1 + P(27:30) * Y(6+P(10)+P(11)) ) ;
     
     % Y_terms(8) = total killing rate applied to V by E [/d]
     % Y_terms(9) = total killing rate applied to V by K [/d]