-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path9dim_model_MAIN.m
101 lines (90 loc) · 2.67 KB
/
9dim_model_MAIN.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
clear;clc;close all
h = 0.01;
t = 0:h:90;
%%
a = 0.6; b = 0.004; c1 = 0.0045; c2 = 0.077; c3 = 0.006; c4 = 0.008;
c5 = 0.045; d1 = 0.5; d2 = 0.8; k1 = 0.6; k2 = 0.5; k3 = 0.6; k4 = 0.65;
%%
e0 = 48; e1 = 0.7; e2 = 2.1; e3 = 1.7; p1 = 0.55; p2 = 0.45; h0 = 270;
h1 = 2.5; h2 = 2; h3 = 10; h4 = 14;
%%
FSH(1) = 142.165;
LH(1) = 28.007;
%%
MsF(1) = 0;
GrF(1) = 0;
PrF(1) = 0;
Sc1(1) = 0;
Sc2(1) = 0;
Lut1(1) = 0;
Lut2(1) = 0;
Lut3(1) = 0;
Lut4(1) = 0;
E2(1) = 50;
P4(1) = 1;
Ih(1) = 290;
% MsF(i) = MsF(i-1) + h * dMsF;
% GrF(i) = GrF(i-1) + h * dGrF;
% PrF(i) = PrF(i-1) + h*dPrF(i);
% Sc1(i) = Scl(i-1) + h*dScl;
% Sc2(i) = Sc2(i-1) + h*dSc2;
% Lut1(i) = Lut1(i-1) + h*dLut1;
% Lut2(i) = Lut2(i-1) + h*dLut2;
% Lut3(i) = Lut3(i-1) + h*dLut3;
% Lut4(i) = Lut4(i-1) + h*dLut4;
%%
for i = 2:length(t)
val = 175*exp(-(t(i)-5)^2/120) + 150*exp(-(t(i)-35)^2/160);
FSH(i) = 250 - (250*(t(i) - 15)^4 )/(1+(t(i)-15)^4) + val;
LH(i) = 380 - 352*(t(i)-15)^4/(1+(t(i)-15)^4);
dMsF = b*FSH(i) + (c1*FSH(i) - c2*LH(i)^a)*MsF(i-1);
MsF(i) = MsF(i-1) + h * dMsF;
dGrF = c2*LH(i)^a*MsF(i-1) + c3*LH(i)^a - c4*LH(i)*GrF(i-1);
GrF(i) = GrF(i-1) + h * dGrF;
dPrF = c4*LH(i)*GrF(i-1) - c5*LH(i)^a*PrF(i-1);
PrF(i) = PrF(i-1) + h*dPrF;
dSc1 = c5*LH(i)^a*PrF(i-1) - d1* Sc1(i-1);
Sc1(i) = Sc1(i-1) + h*dSc1;
dSc2 = d1*Sc1(i-1) - d2* Sc2(i-1);
Sc2(i) = Sc2(i-1) + h*dSc2;
dLut1 = d2*Sc2(i-1) - k1*Lut1(i-1);
Lut1(i) = Lut1(i-1) + h*dLut1;
dLut2 = k1*Lut1(i-1) - k2*Lut2(i-1);
Lut2(i) = Lut2(i-1) + h*dLut2;
dLut3 = k2*Lut2(i-1) - k3*Lut3(i-1);
Lut3(i) = Lut3(i-1) + h*dLut3;
dLut4 = k3*Lut3(i-1) - k4*Lut4(i-1);
Lut4(i) = Lut4(i-1) + h*dLut4;
E2(i) = e0 + e1*GrF(i) + e2*PrF(i) + e3*Lut4(i); %53 nmol/L
P4(i) = p1*Lut3(i) + p2*Lut4(i); %1 nmol/L
Ih(i) = h0 + h1 * PrF(i) + h2 * Lut2(i) + h3 * Lut3(i) + h4 * Lut4(i); %290 U/L
end
figure(1);
subplot(3,1,1); hold on; grid minor
plot(t, E2, 'Linewidth', 2)
xlabel('days')
ylabel('E_2')
title('The amount of Estradiol (E_2) in blood')
subplot(3,1,2); hold on; grid minor
plot(t, P4, 'Linewidth', 2)
xlabel('days')
ylabel('P_4')
title('The amount of Progesterone (P_4) in blood')
subplot(3,1,3); hold on; grid minor
plot(t, Ih, 'Linewidth', 2)
xlabel('days')
ylabel('I_h')
title('The amount of Inhibin (Ih) in blood')
sgtitle('Ovarian Hormones')
figure(2);
subplot(2,1,1); hold on; grid minor
plot(t, FSH, 'Linewidth', 2)
xlabel('days')
ylabel('FSH')
title('The amount of FSH in blood')
subplot(2,1,2); hold on; grid minor
plot(t, LH, 'Linewidth', 2)
xlabel('days')
ylabel('LH')
title('The amount of LH in blood')
sgtitle('Gonadotropin Hormones')