Is there any way to obtain a formula corresponding to this interpolated functions or plots?

I am no expert in this area, but the way I've seen it done is by getting the x,y values from the InterpolatingFunction object, then use InterpolatingPolynomial.

InterpolatingFunction itself uses piecewise InterpolatingPolynomial under the cover. But to obtain just one InterpolatingPolynomial over all the grid points, will result in a very large polynomial of course (order one less that number of grid points), and is not recommended as it will be highly oscillatory. You could pick smaller grid regions at a time, and find the InterpolatingPolynomial over those. But here is what you can try:

ClearAll[x, y, z, gamma, t, a, b, c, d, e];
a = 1; b = 2; c = 3; d = 4; e = 5;
eq1 = x'[t] == -a x[t] - b x[t] y[t] - d x[t] z[t];
eq2 = y'[t] == a x[t] - b y[t] + c z[t] - d x[t] y[t];
eq3 = z'[t] == b y[t] - c z[t] - c x[t] z[t];
eq4 = gamma'[t] == d x[t] y[t] + e x[t] z[t];
sol = First@NDSolve[{eq1, eq2, eq3, eq4, x[0] == 1, y[0] == 0, z[0] == 0, gamma[0] == 0}, 
     {x, y, z, gamma}, {t, 0, 10}];

xSol = x /. sol; (* get formula for the x(t) solution *)
xpoints = Flatten@xSol["Coordinates"];
ypoints = Flatten@xSol["ValuesOnGrid"];
data = Transpose[{xpoints, ypoints}];
p = InterpolatingPolynomial[data, t];
Expand[p]

Here is your formula for the x(t) solution found by NDSolve

DSolve did not give solution:

ClearAll[x, y, z, gamma, t, a, b, c, d, e];
eq1 = x'[t] == -a x[t] - b x[t] y[t] - d x[t] z[t];
eq2 = y'[t] == a x[t] - b y[t] + c z[t] - d x[t] y[t];
eq3 = z'[t] == b y[t] - c z[t] - c x[t] z[t];
eq4 = gamma'[t] == d x[t] y[t] + e x[t] z[t];
DSolve[{eq1, eq2, eq3, eq4}, {x[t], y[t], z[t], gamma[t]}, t]

But NDsolve can solve them. You need to supply numerical values for all the parameters and initial conditions also.

ClearAll[x, y, z, gamma, t, a, b, c, d, e];
a = 1; b = 2; c = 3; d = 4; e = 5;
eq1 = x'[t] == -a x[t] - b x[t] y[t] - d x[t] z[t];
eq2 = y'[t] == a x[t] - b y[t] + c z[t] - d x[t] y[t];
eq3 = z'[t] == b y[t] - c z[t] - c x[t] z[t];
eq4 = gamma'[t] == d x[t] y[t] + e x[t] z[t];
NDSolve[{eq1, eq2, eq3, eq4, x[0] == 1, y[0] == 0, z[0] == 0, 
  gamma[0] == 0}, {x[t], y[t], z[t], gamma[t]}, {t, 0, 10}]