# To reproduce Figure 1 in the associated Physiome paper,
# execute this script from the command line:
#
#   cd [PathToThisFile]
#   [PathToOpenCOR]/pythonshell Figure5.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(B_flow):
    simulation.reset(True)
    simulation.clear_results()
    data.constants()["Apical_concentrations/glucose_m"] = 0.05
    data.constants()["Blood_concentrations/Q_in"] = 1e-17
    data.constants()["phenomonological_constants/n_SGLT"] = 4e7

    data.constants()["Blood_concentrations/Q_in"] = B_flow
    simulation.run()
    ds = simulation.results().data_store()
    glucose_s = ds.voi_and_variables()["Blood_concentrations/glucose_s"].values()[-1] * 1e3
    glucose_i = ds.voi_and_variables()["Cell_concentration/glucose_i"].values()[-1] * 1e3
    J_A_GLUT = ds.voi_and_variables()["A_GLUT2/J_A_GLUT"].values()[-1]
    J_Gl_SGLT = ds.voi_and_variables()["phenomonological_constants/J_Gl_SGLT"].values()[-1]
    J_GLUT = ds.voi_and_variables()["GLUT2/J_GLUT"].values()[-1]

    return (glucose_s, glucose_i, J_A_GLUT, J_Gl_SGLT, J_GLUT)


def run_sim2(n_A_GLUT, n_B_GLUT):
    simulation.reset(True)
    simulation.clear_results()
    data.constants()["Apical_concentrations/glucose_m"] = 0.05
    data.constants()["Blood_concentrations/Q_in"] = 1e-17
    data.constants()["phenomonological_constants/n_SGLT"] = 4e7
    data.constants()["A_GLUT2/n_GLUT"] = n_A_GLUT
    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] * 1e3
    glucose_i = ds.voi_and_variables()["Cell_concentration/glucose_i"].values()[-1] * 1e3
    J_A_GLUT = ds.voi_and_variables()["A_GLUT2/J_A_GLUT"].values()[-1]
    J_Gl_SGLT = ds.voi_and_variables()["phenomonological_constants/J_Gl_SGLT"].values()[-1]
    J_GLUT = ds.voi_and_variables()["GLUT2/J_GLUT"].values()[-1]


    return (glucose_s, glucose_i, J_A_GLUT, J_Gl_SGLT, J_GLUT)


def run_sim3(glucose_int):
    simulation.reset(True)
    simulation.clear_results()
    data.constants()["Apical_concentrations/glucose_m"] = 0.05
    data.constants()["Blood_concentrations/Q_in"] = 1e-17
    data.constants()["phenomonological_constants/n_SGLT"] = 4e7

    # data.constants()["Cell_concentration/L_P"] = 1e-5
    data.constants()["A_GLUT2/n_GLUT"] = 1e8
    data.constants()["GLUT2/n_GLUT"] = 2e8
    data.constants()["Blood_concentrations/glucose_in"] = glucose_int
    simulation.run()
    ds = simulation.results().data_store()
    glucose_s = ds.voi_and_variables()["Blood_concentrations/glucose_s"].values()[-1] * 1e3
    glucose_i = ds.voi_and_variables()["Cell_concentration/glucose_i"].values()[-1] * 1e3
    J_A_GLUT = ds.voi_and_variables()["A_GLUT2/J_A_GLUT"].values()[-1]
    J_Gl_SGLT = ds.voi_and_variables()["phenomonological_constants/J_Gl_SGLT"].values()[-1]
    J_GLUT = ds.voi_and_variables()["GLUT2/J_GLUT"].values()[-1]


    return (glucose_s, glucose_i, J_A_GLUT, J_Gl_SGLT, J_GLUT)




#
#
if __name__ == '__main__':
    # different values for y_shift
    # sglt1 = [1.5e7, 2.25e7, 3e7, 3.75e7, 4.5e7]
    # gl_l = [0.005, 0.01, 0.02, 0.05, 0]
    #
    # X1 = [[1.5e7, 2.25e7, 3e7, 3.75e7, 4.5e7], [0.5e8, 0.75e8, 1e8, 1.25e8, 1.5e8], [1e8, 1.5e8, 2e8, 2.5e8, 3e8],
    #       [1.5e7, 2.25e7, 3e7, 3.75e7, 4.5e7], [0.5e8, 0.75e8, 1e8, 1.25e8, 1.5e8], [1e8, 1.5e8, 2e8, 2.5e8, 3e8],
    #       [1.5e7, 2.25e7, 3e7, 3.75e7, 4.5e7], [0.5e8, 0.75e8, 1e8, 1.25e8, 1.5e8], [1e8, 1.5e8, 2e8, 2.5e8, 3e8],
    #       [1.5e7, 2.25e7, 3e7, 3.75e7, 4.5e7], [0.5e8, 0.75e8, 1e8, 1.25e8, 1.5e8], [1e8, 1.5e8, 2e8, 2.5e8, 3e8],
    #       [1.5e7, 2.25e7, 3e7, 3.75e7, 4.5e7], [0.5e8, 0.75e8, 1e8, 1.25e8, 1.5e8], [1e8, 1.5e8, 2e8, 2.5e8, 3e8]]
    #
    n_A_GLUT = [0, 0.2*2e8, 0.4*2e8, 0.5*2e8]
    n_B_GLUT = [2e8, 0.8*2e8, 0.6*2e8, 0.5*2e8]

    B_flow = [4.5e-18, 1e-17, 2 * 1e-17, 3 * 1e-17]

    n_A_GLUT2 = [0.25*2e8, 0.5*2e8, 0.75*2e8, 1*2e8]
    n_B_GLUT2 = [2e8, 2e8, 2e8, 2e8]

    glucose_int = [0.004, 0.008, 0.012, 0.016]

    y_label = ["Gl$_b$ (mM)", "Gl$_i$ (mM)", "J$_{A_{GLUT2}}$ ($\mu$m/s)", "J$_{SGLT1}$ ($\mu$m/s)",
               "J$_{B_{GLUT2}}$ ($\mu$m/s)"]

    plt.figure(figsize=(15,10))
    plt.subplot(4,5,1)
    Gl_b_results = np.empty(len(B_flow))
    for i in range(len(B_flow)):
        Gl_b = run_sim1(B_flow[i])[0]
        Gl_b_results[i] = Gl_b
    print(Gl_b_results)

    width = 0.5  # the width of the bars
    labels1 = ["0.5Q$_{blood}$", "Q$_{blood}$", "2Q$_{blood}$", "3Q$_{blood}$"]
    plt.bar(labels1, Gl_b_results, width, color='royalblue')
    plt.ylim(0, 16)
    plt.yticks(np.arange(0, 15.1, 5))
    plt.ylabel(y_label[0], fontsize= 12)
    plt.title('A', fontsize= 14)

    plt.subplot(4, 5, 2)
    Gl_i_results = np.empty(len(B_flow))
    for i in range(len(B_flow)):
        Gl_i = run_sim1(B_flow[i])[1]
        Gl_i_results[i] = Gl_i
    print(Gl_i_results)

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

    plt.subplot(4, 5, 3)
    j_AGLUT_results = np.empty(len(B_flow))
    for i in range(len(B_flow)):
        j_AGLUT = run_sim1(B_flow[i])[2]
        j_AGLUT_results[i] = j_AGLUT
    print(j_AGLUT_results)

    width = 0.5  # the width of the bars
    plt.bar(labels1, j_AGLUT_results, width, color='royalblue')
    plt.ylim(0,4e-11)
    plt.yticks(np.arange(0, 3.1e-11, 1e-11))
    plt.ylabel(y_label[2], fontsize=12)
    plt.title('C', fontsize=14)

    plt.subplot(4, 5, 4)
    j_SGLT_results = np.empty(len(B_flow))
    for i in range(len(B_flow)):
        j_SGLT = run_sim1(B_flow[i])[3]
        j_SGLT_results[i] = j_SGLT
    print(j_SGLT_results)

    width = 0.5  # the width of the bars

    plt.bar(labels1, j_SGLT_results, width, color='royalblue')
    plt.ylim(0, 7.4e-11)
    plt.yticks(np.arange(0, 6.1e-11, 2e-11))
    plt.ylabel(y_label[3], fontsize=12)
    plt.title('D', fontsize=14)

    plt.subplot(4, 5, 5)
    j_BGLUT_results = np.empty(len(B_flow))
    for i in range(len(B_flow)):
        j_BGLUT = run_sim1(B_flow[i])[4]
        j_BGLUT_results[i] = j_BGLUT
    print(j_BGLUT_results)

    width = 0.5  # the width of the bars
    plt.bar(labels1, j_BGLUT_results, width, color='royalblue')
    plt.ylim(0, 1.1e-10)
    plt.yticks(np.arange(0, 1.1e-10, 0.2e-10))
    plt.ylabel(y_label[4], fontsize=12)
    plt.title('E', fontsize=14)

    plt.subplot(4, 5, 6)
    Gl_b_results = np.empty(len(n_A_GLUT))
    for i in range(len(n_A_GLUT)):
        Gl_b = run_sim2(n_A_GLUT[i], n_B_GLUT[i])[0]
        Gl_b_results[i] = Gl_b

    width = 0.5  # the width of the bars
    labels2 = ["0:1", "0.2:0.8", "0.4:0.6", "0.5:0.5"]
    plt.bar(labels2, Gl_b_results, width, color='r')
    plt.ylim(0, 12)
    plt.yticks(np.arange(0, 10.1, 2.5))
    plt.ylabel(y_label[0], fontsize=12)
    plt.title('F', fontsize=14)

    plt.subplot(4, 5, 7)
    Gl_i_results = np.empty(len(n_A_GLUT))
    for i in range(len(n_A_GLUT)):
        Gl_i = run_sim2(n_A_GLUT[i], n_B_GLUT[i])[1]
        Gl_i_results[i] = Gl_i

    width = 0.5  # the width of the bars
    labels2 = ["0:1", "0.2:0.8", "0.4:0.6", "0.5:0.5"]
    plt.bar(labels2, Gl_i_results, width, color='r')
    plt.ylim(0, 50)
    plt.yticks(np.arange(0, 40.1, 10))
    plt.ylabel(y_label[1], fontsize=12)
    plt.title('G', fontsize=14)

    plt.subplot(4, 5, 8)
    j_AGLUT_results = np.empty(len(n_A_GLUT))
    for i in range(len(n_A_GLUT)):
        j_AGLUT = run_sim2(n_A_GLUT[i], n_B_GLUT[i])[2]
        j_AGLUT_results[i] = j_AGLUT


    width = 0.5  # the width of the bars

    plt.bar(labels2, j_AGLUT_results, width, color='r')
    plt.ylim(0, 6.2e-12)
    plt.yticks(np.arange(0, 5.1e-12, 1e-12))
    plt.ylabel(y_label[2], fontsize=12)
    plt.title('H', fontsize=14)

    plt.subplot(4, 5, 9)
    j_BGLUT_results = np.empty(len(n_A_GLUT))
    for i in range(len(n_A_GLUT)):
        j_BGLUT = run_sim2(n_A_GLUT[i], n_B_GLUT[i])[3]
        j_BGLUT_results[i] = j_BGLUT

    width = 0.5  # the width of the bars

    plt.bar(labels2, j_BGLUT_results, width, color='r')
    plt.ylim(0, 7.7e-11)
    plt.yticks(np.arange(0, 6.1e-11, 2e-11))
    plt.ylabel(y_label[3], fontsize=12)
    plt.title('I', fontsize=14)

    plt.subplot(4, 5, 10)
    Gl_b_results = np.empty(len(n_A_GLUT))
    for i in range(len(n_A_GLUT)):
        Gl_b = run_sim2(n_A_GLUT[i], n_B_GLUT[i])[4]
        Gl_b_results[i] = Gl_b

    width = 0.5  # the width of the bars

    plt.bar(labels2, Gl_b_results, width, color='r')
    plt.ylim(0, 8e-11)
    plt.yticks(np.arange(0, 6.1e-11, 2e-11))
    plt.ylabel(y_label[4], fontsize=12)
    plt.title('I', fontsize=14)

    plt.subplot(4, 5, 11)
    Gl_b_results = np.empty(len(n_A_GLUT2))
    for i in range(len(n_A_GLUT2)):
        Gl_b = run_sim2(n_A_GLUT2[i], n_B_GLUT2[i])[0]
        Gl_b_results[i] = Gl_b

    width = 0.5  # the width of the bars
    labels3 = ["0.25:1", "0.5:1", "0.75:1", "1:1"]
    plt.bar(labels2, Gl_b_results, width, color='green')
    plt.ylim(0, 13)
    plt.yticks(np.arange(0, 12.6, 2.5))
    plt.ylabel(y_label[0], fontsize=12)
    plt.title('K', fontsize=14)

    plt.subplot(4, 5, 12)
    Gl_i_results = np.empty(len(n_A_GLUT2))
    for i in range(len(n_A_GLUT2)):
        Gl_i = run_sim2(n_A_GLUT2[i], n_B_GLUT2[i])[1]
        Gl_i_results[i] = Gl_i

    width = 0.5  # the width of the bars
    plt.bar(labels3, Gl_i_results, width, color='green')
    plt.ylim(0, 40)
    plt.yticks(np.arange(0, 30.1, 10))
    plt.ylabel(y_label[1], fontsize=12)
    plt.title('L', fontsize=14)

    plt.subplot(4, 5, 13)
    j_AGLUT_results = np.empty(len(n_A_GLUT2))
    for i in range(len(n_A_GLUT2)):
        j_AGLUT = run_sim2(n_A_GLUT2[i], n_B_GLUT2[i])[2]
        j_AGLUT_results[i] = j_AGLUT

    width = 0.5  # the width of the bars

    plt.bar(labels3, j_AGLUT_results, width, color='green')
    plt.ylim(0, 2.3e-11)
    plt.yticks(np.arange(0, 2.1e-11, 0.5e-11))
    plt.ylabel(y_label[2], fontsize=12)
    plt.title('M', fontsize=14)

    plt.subplot(4, 5, 14)
    j_BGLUT_results = np.empty(len(n_A_GLUT2))
    for i in range(len(n_A_GLUT2)):
        j_BGLUT = run_sim2(n_A_GLUT2[i], n_B_GLUT2[i])[3]
        j_BGLUT_results[i] = j_BGLUT

    width = 0.5  # the width of the bars

    plt.bar(labels3, j_BGLUT_results, width, color='green')
    plt.ylim(0, 7.5e-11)
    plt.yticks(np.arange(0, 6.1e-11, 2e-11))
    plt.ylabel(y_label[3], fontsize=12)
    plt.title('N', fontsize=14)

    plt.subplot(4, 5, 15)
    Gl_b_results = np.empty(len(n_A_GLUT2))
    for i in range(len(n_A_GLUT2)):
        Gl_b = run_sim2(n_A_GLUT2[i], n_B_GLUT2[i])[4]
        Gl_b_results[i] = Gl_b

    width = 0.5  # the width of the bars

    plt.bar(labels3, Gl_b_results, width, color='green')
    plt.ylim(0, 9e-11)
    plt.yticks(np.arange(0, 8.1e-11, 2e-11))
    plt.ylabel(y_label[4], fontsize=12)
    plt.title('O', fontsize=14)

    plt.subplot(4, 5, 16)
    Gl_b_results = np.empty(len(glucose_int))
    for i in range(len(glucose_int)):
        Gl_b = run_sim3(glucose_int[i])[0]
        Gl_b_results[i] = Gl_b
    print(Gl_b_results)

    width = 0.5  # the width of the bars
    labels4 = ["4", "8", "12", "16"]
    plt.bar(labels4, Gl_b_results, width, color='orange')
    plt.ylim(0, 21)
    plt.yticks(np.arange(0, 20.1, 5))
    plt.ylabel(y_label[0], fontsize=12)
    plt.title('P', fontsize=14)

    plt.subplot(4, 5, 17)

    Gl_i_results = np.empty(len(glucose_int))
    for i in range(len(glucose_int)):
        Gl_i = run_sim3(glucose_int[i])[1]
        Gl_i_results[i] = Gl_i

    width = 0.5  # the width of the bars
    plt.bar(labels4, Gl_i_results, width, color='orange')
    plt.ylim(0, 52)
    plt.yticks(np.arange(0, 40.1, 10))
    plt.ylabel(y_label[1], fontsize=12)
    plt.title('Q', fontsize=14)

    plt.subplot(4, 5, 18)
    j_AGLUT_results = np.empty(len(glucose_int))
    for i in range(len(glucose_int)):
        j_AGLUT = run_sim3(glucose_int[i])[2]
        j_AGLUT_results[i] = j_AGLUT

    width = 0.5  # the width of the bars

    plt.bar(labels4, j_AGLUT_results, width, color='orange')
    plt.ylim(0, 1.8e-11)
    plt.yticks(np.arange(0, 1.51e-11, 0.5e-11))
    plt.ylabel(y_label[2], fontsize=12)
    plt.title('R', fontsize=14)

    plt.subplot(4, 5, 19)
    j_BGLUT_results = np.empty(len(glucose_int))
    for i in range(len(glucose_int)):
        j_BGLUT = run_sim3(glucose_int[i])[3]
        j_BGLUT_results[i] = j_BGLUT

    width = 0.5  # the width of the bars

    plt.bar(labels4, j_BGLUT_results, width, color='orange')
    plt.ylim(0, 7.5e-11)
    plt.yticks(np.arange(0, 6.1e-11, 2e-11))
    plt.ylabel(y_label[3], fontsize=12)
    plt.title('S', fontsize=14)

    plt.subplot(4, 5, 20)
    Gl_b_results = np.empty(len(glucose_int))
    for i in range(len(glucose_int)):
        Gl_b = run_sim3(glucose_int[i])[4]
        Gl_b_results[i] = Gl_b

    width = 0.5  # the width of the bars

    plt.bar(labels4, Gl_b_results, width, color='orange')
    plt.ylim(0, 9e-11)
    plt.yticks(np.arange(0, 8.1e-11, 2e-11))
    plt.ylabel(y_label[4], fontsize=12)
    plt.title('T', fontsize=14)
    plt.tight_layout()


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