(*NOTE:All concentrations are in M.*)(*Constants.*)(*Note:R is in \
units of J/(mol*K).T is in units of Kelvin.*)
R = 8.314;
T = 310.15;

(*Kd values for protein-oligo interactions and delta G naught values \
for oligo-oligo interactions.Values of Kd are in mol/L.Delta G naught \
values are in joules/mol*)
KdA1P = 0;
KdA2P = 0;
GnaughtA1C = -46860.8;
GnaughtA2C = -37237.6; (*or \
29706.4 for 2C or 22593.6 for 1C connectors*)

(*Values used for mass conservation equations.Equal to the initial \
concentrations of each molecule.*)
A1total = 1;
A2total = 1;
Ctotal = 0.5;
Ptotal = 0.1;

(*Mass balance equations*)
massBalance1 =   A1total == A1 + A1P + A1A2P + A1C + A1A2C + A1PC + A1A2PC;
massBalance2 =   A2total == A2 + A2P + A1A2P + A2C + A1A2C + A2PC + A1A2PC;
massBalance3 = Ptotal == P + A1P + A2P + A1A2P + A1PC + A2PC + A1A2PC;
massBalance4 =   Ctotal == Connector + A1C + A2C + A1A2C + A1PC + A2PC + A1A2PC;
massBalanceEquations = {massBalance1, massBalance2, massBalance3,  massBalance4};

(*Concentrations cannot be negative*)
nonNegativeConstraint = {A1 >= 0, A2 >= 0, Connector >= 0, P >= 0, A1P >= 0, A2P >= 0, A1A2P >= 0, A1C >= 0, A2C >= 0, A1A2C >= 0, A1PC >= 0, A2PC >= 0, A1A2PC >= 0};

(*this constraint is not strictly necessary but may improve speed of \
the minimization function by limiting the range of possible values*)

maxValueConstraint = {A1 <= A1total + .001, A2 <= A2total + .001, Connector <= Ctotal + .001, P <= Ptotal + .001, A1P <= Ptotal + .001, A2P <= Ptotal + .001, A1A2P <= Ptotal + .001, A1PC <= Ptotal + .001, A2PC <= Ptotal + .001, A1A2PC <= Ptotal + .001, A1C <= Ctotal + .001, A2C <= Ctotal + .001, A1A2C <= Ctotal + .001};

(*deltaG=-R*T*Log[1/Kd]+R*T*Log[()/()]*)
deltaGA1P = Abs[-R*T*Log[1/KdA1P] + R*T*Log[(A1P)/(A1*P)]];
deltaGA2P = Abs[-R*T*Log[1/KdA2P] + R*T*Log[(A2P)/(A2*P)]];
deltaGA1A2P = Abs[-R*T*Log[1/KdA1P] - R*T*Log[1/KdA2P] + R*T*Log[(A1A2P)/(A1P*A2)] + R*T*Log[(A1A2P)/(A2P*A1)]];
deltaGA1C = Abs[GnaughtA1C + R*T*Log[(A1C)/(A1*Connector)]];
deltaGA2C = Abs[GnaughtA2C + R*T*Log[(A2C)/(A2*Connector)]];
deltaGA1A2C = Abs[GnaughtA1C + GnaughtA2C + R*T*Log[(A1A2C)/(A1C*A2)] + R*T*Log[(A1A2C)/(A2C*A1)]];
deltaGA1PC = Abs[GnaughtA1C - R*T*Log[1/KdA1P] + R*T*Log[(A1PC)/(A1P*Connector)] + R*T*Log[(A1PC)/(A1C*P)]];
deltaGA2PC = Abs[GnaughtA2C - R*T*Log[1/KdA2P] +  R*T*Log[(A2PC)/(A2P*Connector)] + R*T*Log[(A2PC)/(A2C*P)]];
deltaGA1A2PC = Abs[(GnaughtA1C + GnaughtA2C) + R*T*Log[(A1A2PC)/(A1A2P*Connector)] + (GnaughtA1C - R*T*Log[1/KdA1P]) + R*T*Log[(A1A2PC)/(A2PC*A1)] + (GnaughtA2C - R*T*Log[1/KdA2P]) + R*T*Log[(A1A2PC)/(A1PC*A2)] + (-R*T*Log[1/KdA1P] - R*T*Log[1/KdA2P]) + R*T*Log[(A1A2PC)/(A1A2C*P)] + (GnaughtA1C - R*T*Log[1/KdA2P]) + R*T*Log[(A1A2PC)/(A1P*A2C)] + (GnaughtA2C - R*T*Log[1/KdA1P]) + R*T*Log[(A1A2PC)/(A2P*A1C)]];

(*Equation for the total delta G of the full reaction.The value of \
this will be at a minimum when the reaction reaches equilibrium.*)
sumOfAllDeltaG = deltaGA1C + deltaGA2C + deltaGA1A2C + deltaGA1P + deltaGA2P + 
   deltaGA1A2P + deltaGA1PC + deltaGA2PC + deltaGA1A2PC;

(*Uses the minimize function to find the concentrations of all \
molecules and complexes at equilibrium given the constraints.*)
NMinimize[{sumOfAllDeltaG, {massBalanceEquations, nonNegativeConstraint (*, maxValueConstraint*)}}, {A1, A2, Connector, P, A1P, A2P, A1A2P, A1C, A2C, A1A2C, A1PC, A2PC, A1A2PC}, WorkingPrecision -> 1000, PrecisionGoal -> 1000, AccuracyGoal -> 1000]
Clear["Global`*"]