diff --git a/Runge_Kutta_solution.py b/Runge_Kutta_solution.py new file mode 100644 index 0000000000000000000000000000000000000000..e1c0a4ca6ab99bfe30add79719ad39d1d2347fa9 --- /dev/null +++ b/Runge_Kutta_solution.py @@ -0,0 +1,76 @@ +import numpy as np +import matplotlib.pyplot as plt + +plt.rcParams['font.size'] = 14 +plt.rcParams['font.family'] = 'Times New Roman' + +def dydt(t, Y, Kp, Km, X0, A0): + return Kp * (X0 - Y) * (A0 - Y) - Km * Y + +def runge_kutta(h, t_max, Kp, Km, X0, A0): + t_values = np.arange(0, t_max + h, h) + y_values = np.zeros(len(t_values)) + y = 0 + + for i in range(1, len(t_values)): + t = t_values[i-1] + k1 = h * dydt(t, y, Kp, Km, X0, A0) + k2 = h * dydt(t + 0.5*h, y + 0.5*k1, Kp, Km, X0, A0) + k3 = h * dydt(t + 0.5*h, y + 0.5*k2, Kp, Km, X0, A0) + k4 = h * dydt(t + h, y + k3, Kp, Km, X0, A0) + + y += (k1 + 2*k2 + 2*k3 + k4) / 6 + y_values[i] = y + + return t_values, y_values + + +#Set the parameters. +#Association constant k+. +Kp = @Kp +#Dissociation constant k-. +Km = @Km +#Initial free binding sites X0 (receptors or antigens). +X0 = @Xo +#Initial AB concentration A0. +A0 = @A0 +#Step width of Runge-Kutta algorythm. Can be set to 0.1 for rough estimation. +h = @h0 +#Measuring time t_max. +t_max = @t_max + +#Insert time points. +t_data = @t_data #@t_data=np.array([t1,t2,t3,...]) +#Insert MFI data. +y_data = @y_data #@y_data=np.array([y1,y2,y3,...]) +#Insert y standard deviation. +stdev_data = @stdev_data #@stdev_data=np.array([stdev1,stdev2,stdev3,...]) + + +t_values, y_values = runge_kutta(h, t_max, Kp, Km, X0, A0) + +plt.plot(t_data, y_data, 'o', color='green', label='Label') + + +#Plot the results. +plt.plot(t_values, y_values, color='#8A2BE2', label='Runge-Kutta Solution') +lower_bound = y_data - stdev_data +upper_bound = y_data + stdev_data +plt.fill_between(t_data, lower_bound, upper_bound, color='lightgreen', alpha=0.3, label='SD') +plt.xlabel('Time in s') +plt.ylabel('MFI') +plt.legend(prop={'weight': 'bold'}, loc='lower right') +plt.title('Title', fontweight='bold') + +#Plot limits. +plt.xlim(-5, 3600) +plt.ylim(0,150) + +ax = plt.gca() +ax.spines['top'].set_visible(False) +ax.spines['right'].set_visible(False) + + +#plt.savefig('Figure.png', bbox_inches='tight') + +plt.show() \ No newline at end of file