adherence_project/RK4.cpp at master · JDChallenger/adherence_project · GitHub

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

#include "parameters.h"
#include <vector>
#include <cmath>
#include <cstdlib>

using namespace std;

//////////////////////////////////////////////////
/* Define RK4 functions */
//////////////////////////////////////////////////
double f1AM(double x,double y, double z , params* theta)//(x,y,z)=(gut,central,metabolite)
{
double temp;
temp = -(theta->kaAM) * x;
return temp;
}

double f2AM(double x,double y, double z , params* theta)//(x,y,z)=(gut,central,metabolite)
{
double temp;
temp = (theta->kaAM) * x -((theta->k23AM) + ((theta->CLAM) / (theta->VCAM)))* y;
return temp;
}

double f3AM(double x,double y, double z , params* theta)//(x,y,z)=(gut,central,metabolite)
{
double temp;
temp = (theta->k23AM) * y - ((theta->CLmetAM) / (theta->VMAM))*z;
return temp;
}

double f1L(double x,double y, double z , params* theta)//(x,y,z)=(gut,central,metabolite)
{
double temp;
temp = -(theta->kaL) * x;
return temp;
}

double f2L(double x,double y, double z , params* theta)//(x,y,z)=(gut,central,metabolite)
{
double temp;
temp = (theta->kaL) * x -((theta->k23L) + ((theta->CLL) / (theta->VCL)))* y;
return temp;
}

double f3L(double x,double y, double z, params* theta)//(x,y,z)=(gut,central,metabolite)
{
double temp;
temp = (theta->k23L) * y - ((theta->CLmetL) / (theta->VML))*z;
return temp;
}

//Implementation of the analytical solution for the Artemether PK equation
vector<double> ART(double Time, double Dose, params* theta){
	vector<double> temp(3);

	double c1 = - (Dose * theta->k23AM * theta->kaAM * theta->VCAM * theta->VMAM * theta->VMAM)/
	((-theta->CLmetAM + theta->kaAM*theta->VMAM)*
		(-theta->CLmetAM * theta->VCAM + theta->CLAM * theta->VMAM + theta->k23AM * theta->VCAM * theta->VMAM));
	double c2 = (Dose * theta->k23AM * theta->kaAM * theta->VCAM * theta->VMAM)/
	((-theta->CLmetAM + theta->kaAM*theta->VMAM)*
		(theta->CLAM + theta->k23AM * theta->VCAM - theta->kaAM * theta->VCAM));
	double c3 = -(Dose * theta->k23AM * theta->kaAM * theta->VCAM * theta->VCAM * theta->VMAM)/
	((-theta->CLmetAM * theta->VCAM + theta->CLAM * theta->VMAM + theta->k23AM * theta->VCAM * theta->VMAM)*
		(theta->CLAM + theta->k23AM * theta->VCAM - theta->kaAM * theta->VCAM));

	double eig1 = -theta->CLmetAM / theta->VMAM;
	double eig2 = -theta->kaAM;
	double eig3 = -(theta->CLAM + theta->k23AM * theta->VCAM)/(theta->VCAM);

	double eigV1[3] = {0,0,1};
	double eigV2[3] =
		{((theta->CLAM + theta->k23AM * theta->VCAM - theta->kaAM * theta->VCAM)*
			(-theta->CLmetAM + theta->kaAM * theta->VMAM))/(theta->k23AM * theta->kaAM * theta->VCAM * theta->VMAM),
		-(theta->CLmetAM - theta->kaAM * theta->VMAM)/(theta->k23AM * theta->VMAM),1};
	double eigV3[3] = {0,
		-(theta->CLmetAM * theta->VCAM - theta->CLAM * theta->VMAM - theta->k23AM * theta->VCAM * theta->VMAM)/
		(theta->k23AM * theta->VCAM * theta->VMAM),1};

	//Eigenvector #3 always has a factor of -1 which we need to remove
	temp[0] = c1 * exp(eig1*Time) * eigV1[0] + c2 * exp(eig2*Time) * eigV2[0] + c3 * exp(eig3*Time) * eigV3[0];
	temp[1] = c1 * exp(eig1*Time) * eigV1[1] + c2 * exp(eig2*Time) * eigV2[1] + c3 * exp(eig3*Time) * eigV3[1];
	temp[2] = abs(c1 * exp(eig1*Time) * eigV1[2] + c2 * exp(eig2*Time) * eigV2[2] + c3 * exp(eig3*Time) * eigV3[2]);

	return temp;
}