dx1/dt = (-ka*n*x1 - ka*kc*a1*x1 + ka*kc*a1*x2 - ka*a1*x1^2 + ka*a1*x1*x2 - kc*a1*x1*x3 + kc*a1*x2*x3) / (ka*kc + ka*x1 + kc*x3); 
dx2/dt = a2*x1 - a2*x2;
dx3/dt = (ka*kc*b1*x4 - ka*kc*b1*x3 + ka*b1*x1*x4 - ka*b1*x1*x3 - n*kc*x3 + kc*b1*x4*x3 - kc*b1*x3^2) / (ka*kc + ka*x1 + kc*x3);  
dx4/dt = -b2*x4 + b2*x3;
y1 = x1
