This model has two state variables, X1 and X2, with dynamics
dX1 = \[Mu]1[x1,x2] dt + \[Sigma]11[x1,x2] dW1 + \[Sigma]12[x1,x2] dW2
dX2 = \[Mu]2[x1,x2] dt + \[Sigma]21[x1,x2] dW1 + \[Sigma]22[x1,x2] dW2
The two Brownian motions are independent; any correlation comes through
the off-diagonal terms in the \[Sigma] matrix.
Below is an example, with parameters a1, b1, a2, b2, g11, g22, g21, r:
\[Mu]1[x1_,x2_] = a1 + b1*x1;
\[Mu]2[x1_,x2_] = a2 + b2*x2;
\[Sigma]11[x1_,x2_] = g11*r*Sqrt[x1];
\[Sigma]22[x1_,x2_] = g22*Sqrt[1-r^2]*Sqrt[x1];
\[Sigma]12[x1_,x2_] = 0;
\[Sigma]21[x1_,x2_] = g21;
Edit these expressions to match your model, and email to me.
Remark: do not use the letters c and k in your parameter names as they
conflict with other symbols in the code.