#WRITE THE MODEL STRING HERE AND LOAD IT TO TELLURIUM
model_stringC4='''
  # Identify variables and specify initial values:
  AcP = 0; HOAc = 30; ACoA = 0; Pta = 0; Acs = 0; LacI = 0;
 
  # Reactions with rate laws:
  R1: -> ACoA ; k1  
  R2: ACoA -> ; k2*ACoA
  R3: ACoA -> AcP; k3*Pta*ACoA/(Kpta + ACoA)
  R4: AcP -> HOAc; kf*AcP
  R5: HOAc -> AcP; kr*HOAc
  R6: HOAc -> ACoA; k4*Acs*HOAc/(Kacs+HOAc)
  R7: HOAc ->; k5*HOAc
  R8: -> Acs; k6*(AcP^n1)/(Kacp^n1 + AcP^n1)
  R9: -> LacI; k7*(AcP^n2)/(Kacp^n2 + AcP^n2)
  R10: -> Pta; Vmax*(1-(LacI^n3/(LacI^n3 + Krpta^n3)))
  R11: Pta -> ; k9*Pta
  R12: Acs -> ; k9*Acs
  R13: LacI ->; k9*LacI
  R14: -> HOAc; k10
  # Kinetic parameters:
 k1 = 0.5; k2 = 10; k3 = 80000; kf = 1; kr = 1; k4 = 800; k5 = 0.005; k6 = 0.002; k7 = 0.0001; Vmax = 0.002; k9 = 0.06; k10 = 0.005 ;n1 = 2; n2 = 2; n3 = 2;
 Kpta = 0.06; Kacs = 0.1; Kacp = 10; Krpta = 0.000001
  '''
 
modelC4 = te.loada(model_stringC4)

#SIMULATE THE MODEL HERE
modelC4.resetAll()
resultC4 = modelC4.simulate(start=0, end = 3000, points=10000)

#REPLACE NAMES OF THE VARIABLES USED HERE TO ONES YOU USED IN THE SIMULATION

plt.figure(1, (8, 8))

ax = plt.subplot(2, 2, 1)
plt.plot(resultC4['time'], resultC4['[AcP]'],label='AcP')
plt.plot(resultC4['time'], resultC4['[HOAc]'], label='HOAc')

plt.legend()
plt.xlabel('Time, min')
plt.ylabel('Concentration, mM')

ax = plt.subplot(2, 2, 2)
plt.plot(resultC4['time'], resultC4['[ACoA]'],label='ACoA',color="red")
plt.legend()
plt.xlabel('Time, min')
plt.ylabel('Concentration, mM')

ax = plt.subplot(2, 2, 3)
plt.plot(resultC4['time'], resultC4['[Pta]'],label='Pta',color="purple")

plt.legend()
plt.xlabel('Time, min')
plt.ylabel('Concentration, mM')

ax = plt.subplot(2, 2, 4)
plt.plot(resultC4['time'], resultC4['[Acs]'],label='Acs',color="violet")
plt.plot(resultC4['time'], resultC4['[LacI]'], label='LacI',color="darkred")

plt.legend()
plt.xlabel('Time, min')
plt.ylabel('Concentration, mM')

plt.tight_layout()