# To reproduce Figure 1 in the associated Physiome paper,
# execute this script from the command line:
#
#   cd [PathToThisFile]
#   [PathToOpenCOR]/pythonshell Figure6.py

import matplotlib.pyplot as plt
import opencor as opencor
import numpy as np

simulation = opencor.open_simulation("model.sedml")
data = simulation.data()
data.set_ending_point(3000)
data.set_point_interval(10)




def run_sim1(n_A_GLUT, n_B_GLUT, n_SGLT1):
    simulation.reset(True)
    simulation.clear_results()
    data.constants()["Apical_concentrations/glucose_m"] = 0.05
    data.constants()["Cell_concentration/L_A"] = 5e-5
    data.constants()["Cell_concentration/L_B"] = 5e-5
    data.constants()["Blood_concentrations/v_B"] = 1e-17
    data.constants()["Blood_concentrations/glucose_in"] = 0.005
    data.constants()["Blood_concentrations/Q_in"] = 9e-18
    data.constants()["A_GLUT2/n_GLUT"] = n_A_GLUT
    data.constants()["phenomonological_constants/n_SGLT"] = n_SGLT1
    data.constants()["GLUT2/n_GLUT"] = n_B_GLUT
    simulation.run()
    ds = simulation.results().data_store()
    glucose_s = ds.voi_and_variables()["Blood_concentrations/glucose_s"].values()[-1]
    glucose_i = ds.voi_and_variables()["Cell_concentration/glucose_i"].values()[-1]
    J_GLUT = ds.voi_and_variables()["GLUT2/J_GLUT"].values()[-1]
    J_Na_NaK = ds.voi_and_variables()["NaK_channel/J_Na_NaK"].values()[-1]

    return (glucose_s, glucose_i, J_GLUT, J_Na_NaK)


#
if __name__ == '__main__':

    #
    n_A_GLUT = [0, 0.25*2e8, 0.5*2e8, 0.25*2e8, 1e5]
    n_B_GLUT = [2e8, 0.75*2e8, 0.5*2e8, 0.75*2e8, 2e8]
    n_SGLT1 = [4e7, 6e7, 4e7, 8e7, 8e7]

    y_label = ["Blood Glucose", "Cell Glucose", "GLUT2 Flux"]

    plt.figure(figsize=(15,10))
    plt.subplot(1,3,1)
    Gl_b_results = np.empty(len(n_A_GLUT))
    for i in range(len(n_A_GLUT)):
        Gl_b = run_sim1(n_A_GLUT[i], n_B_GLUT[i], n_SGLT1[i])[0]
        Gl_b_results[i] = Gl_b
    print(Gl_b_results)

    width = 0.5  # the width of the bars
    labels1 = ["A", "B", "C", "D", "E"]
    plt.bar(labels1, Gl_b_results, width, color='r')
    plt.ylim(0, 0.014)
    plt.yticks(np.arange(0, 0.0121, 0.002))
    plt.ylabel(y_label[0], fontsize= 12)
    plt.title('A', fontsize= 14)

    plt.subplot(1, 3, 2)
    Gl_i_results = np.empty(len(n_A_GLUT))
    for i in range(len(n_A_GLUT)):
        Gl_i = run_sim1(n_A_GLUT[i], n_B_GLUT[i], n_SGLT1[i])[1]
        Gl_i_results[i] = Gl_i

    width = 0.5  # the width of the bars
    plt.bar(labels1, Gl_i_results, width, color='b')
    plt.ylim(0, 0.063)
    plt.yticks(np.arange(0, 0.061, 0.02))
    plt.ylabel(y_label[1], fontsize=12)
    plt.title('B', fontsize=14)

    plt.subplot(1, 3, 3)
    J_GLUT_results = np.empty(len(n_A_GLUT))
    for i in range(len(n_A_GLUT)):
        J_GLUT = run_sim1(n_A_GLUT[i], n_B_GLUT[i], n_SGLT1[i])[2]
        J_GLUT_results[i] = J_GLUT

    width = 0.5  # the width of the bars
    plt.bar(labels1, J_GLUT_results, width, color='g')
    plt.ylim(0, 10e-11)
    plt.yticks(np.arange(0, 8.1e-11, 2e-11))
    plt.ylabel(y_label[2], fontsize=12)
    plt.title('C', fontsize=14)

    plt.savefig('Figure06.png')
    plt.show()