function xprime = prostateON(t, x)
% System of ODEs with prostate cancer dynamics
% Androgen suppression is ON

% Constants (values from Table 1 in [Liu2015])
% Patient 1
alpha1 = 0.0204;
alpha2 = 0.0242;
beta1 = 0.0201;
beta2 = 0.0168;
k1 = 10.0*10^-9;
k2 = 1.0;
k3 = 10.0*10^-9;
k4 = 2.0;
m1 = 0.00005;
x30 = 12.0;
tau = 12.5;
lambda1 = 0.03;
%lambda1 = 0.01;
mu1 = 0.05;
mu3 = 0.02;
d = 1.0;

% Patient 15
% alpha1 = 0.0213;
% alpha2 = 0.0242;
% beta1 = 0.01;
% beta2 = 0.0168;
% k1 = 7*10^-9;
% k2 = 1.0;
% k3 = 7.4*10^-9;
% k4 = 2.0;
% m1 = 0.00005;
% x30 = 8.0;
% tau = 12.5;
% lambda1 = 0.03;
% mu1 = 0.02;
% mu3 = 0.02;
% d = 1.0;

% Patient 26
% alpha1 = 0.0197;
% alpha2 = 0.0242;
% beta1 = 0.0175;
% beta2 = 0.0168;
% k1 = 10*10^-9;
% k2 = 1.0;
% k3 = 10*10^-9;
% k4 = 2.0;
% m1 = 0.00005;
% x30 = 12.0;
% tau = 12.5;
% lambda1 = 0.01;
% mu1 = 0.03;
% mu3 = 0.02;
% d = 1.0;

% Patient 11
% alpha1 = 0.0204;
% alpha2 = 0.0242;
% beta1 = 0.02;
% beta2 = 0.0158;
% k1 = 7*10^-9;
% k2 = 1.0;
% k3 = 7*10^-9;
% k4 = 2.0;
% m1 = 0.00005;
% x30 = 9.0;
% tau = 12.5;
% lambda1 = 0.0121;
% mu1 = 0.06;
% mu3 = 0.02;
% d = 1.0;

% theta3 = 0.0213; % HSC proliferation constant (alpha1)
% theta4 = 0.0242; % CRC proliferation constant (alpha2)
% theta5 = 0.01; % HSC apoptosis constant (beta1)
% theta6 = 0.0168; % CRC apoptosis constant (beta2)
% 
% alpha1 = theta3;
% alpha2 = theta4;
% beta1 = theta5;
% beta2 = theta6;

% Without white noise
xprime = [(alpha1/(1+exp(-(x(3)-k1)*k2)))*x(1) - (beta1/(1+exp(-(x(3)-k3)*k4)))*x(1)...
    - m1*(1 - (x(3)/x30))*x(1) - lambda1*x(1) + mu1;...
    (alpha2*(1- (d*x(3)/x30)))*x(2) - beta2*x(2) + m1*(1 - (x(3)/x30))*x(1);...
    - (1/tau)*x(3) + mu3];

% % Parameters of Gaussian white noise
% mymean = 0.0;
% mystdeviation = 0.001; %sqrt(0.001);
% 
% % With white noise in x1 and x2
% xprime = [(alpha1/(1+exp(-(x(3)-k1)*k2)))*x(1) - (beta1/(1+exp(-(x(3)-k3)*k4)))*x(1)...
%     - m1*(1 - (x(3)/x30))*x(1) - lambda1*x(1) + mu1 + (mystdeviation*randn(1)+mymean);...
%     (alpha2*(1- (d*x(3)/x30)))*x(2) - beta2*x(2) + m1*(1 - (x(3)/x30))*x(1) + (mystdeviation*randn(1)+mymean);...
%     - (1/tau)*x(3) + mu3];

% With white noise in x1, x2 and x3
% xprime = [(alpha1/(1+exp(-(x(3)-k1)*k2)))*x(1) - (beta1/(1+exp(-(x(3)-k3)*k4)))*x(1)...
%     - m1*(1 - (x(3)/x30))*x(1) - lambda1*x(1) + mu1 + (mystdeviation*randn(1)+mymean);...
%     (alpha2*(1- (d*x(3)/x30)))*x(2) - beta2*x(2) + m1*(1 - (x(3)/x30))*x(1) + (mystdeviation*randn(1)+mymean);...
%     - (1/tau)*x(3) + mu3 + (mystdeviation*randn(1)+mymean)];

end