%% Change in free energy
% assume the temperature is 25 degree
koff = 2.7e-2; %1/s
kon = 1.9e4 % 1/Ms
Kd = kon/koff
Ka = 1/Kd
% Keq = Kd = 1/Ka
R = 8.314472; % J/K*mol universal gas constant
T = 25+273; % K
dG = log(Kd)*R*T % change in free energy
%% number of free and occupied binding site
S = 10e-6; %amount of streptavidin (μM)
B = 10e-6; %amount of biotin (μM)
S2 = S+B
S3 = S+2*B
S4 = S+3*B
S5 = S+4*B
syms S(t) S2(t) S3(t) S4(t) S5(t)
ode1 = diff(S,t) == (-4*kon*B*S)+koff*S2;
ode2 = diff(S2,t) == -(3*kon*B+koff)*S2+4*kon*B*S+2*koff*S3;
ode3 = diff(S3,t) == -(2*kon*B+2*koff)*S3+3*kon*B*S2+3*koff*S4;
ode4 = diff(S4,t) == -(kon*B+3*koff)*S4+2*kon*B*S3+4*koff*S5;
ode5 = diff(S5,t) == -4*koff*S5+kon*B*S4;
odes = [ode1; ode2; ode3; ode4; ode5]
Aodes = dsolve(odes)
Here's the full version of the code.