lst = Table[X0 = 600;
  Y0 = 300;  [Mu]x = 0.006; [Mu]y = 0.006; [Mu]2 =
   0.00001; [Lambda] = 1; [Sigma]x = 1;  [Sigma]y = 1;  [Gamma]x =
   0.1; [Gamma]y = 0.2; [Epsilon]y = 10; [Phi] = 10; [Alpha] =
   0.008;
  gx[t_, [Phi]_, [Lambda]_, [Epsilon]x_] :=
   If[t < [Epsilon]x, 0,
    PDF[GammaDistribution[[Phi], 1/[Lambda]], t - [Epsilon]x]];
  gy[t_, [Phi]_, [Lambda]_, [Epsilon]y_] :=
   If[t < [Epsilon]y, 0,
    PDF[GammaDistribution[[Phi], 1/[Lambda]], t - [Epsilon]y]];
  soln[[Epsilon]x_] =
   sol := NDSolve[{
      !(
*SubscriptBox[([PartialD]), (t)](x1[t])) ==
       x3[t] gx[t, [Phi], [Lambda], [Epsilon]x] -
        x1[t] ([Gamma]x + [Mu]x),
      !(
*SubscriptBox[([PartialD]), (t)](x2[
         t])) == [Gamma]x x1[t] - [Mu]2 x2[t],
      !(
*SubscriptBox[([PartialD]), (t)](y1[t])) ==
       y3[t] gy[t, [Phi], [Lambda], [Epsilon]y] -
        y1[t] ([Gamma]y + [Mu]y),
      !(
*SubscriptBox[([PartialD]), (t)](y2[
         t])) == [Gamma]y y1[t] - [Mu]2 y2[t],
      !(
*SubscriptBox[([PartialD]), (t)](x3[t])) == 0,
      !(
*SubscriptBox[([PartialD]), (t)](y3[t])) == 0,
      x1[0] == 0, x2[0] == 0, y1[0] == 0, y2[0] == 0, x3[0] == X0,
      y3[0] == Y0},
     {x1, x2, x3, y1, y2, y3}, {t, 0, 120}];
  Xtotal = {X0}; Ytotal = {Y0};
  For[y = 0, y < 30, y++,
   loop = sol;
   X0 = [Sigma]x y2[120] /. loop[[1, 4]];
   Y0 = [Sigma]y x2[120]/(1 + [Alpha] x2[120]) /. loop[[1, 3]];
   Xtotal = Append[Xtotal, X0];
   Ytotal = Append[Ytotal, Y0];
   ];, {[Epsilon]x, 0, 10, 5}]
ListPlot[Xtotal]
ListPlot[Ytotal]