Manipulate an integral with variable in limit of the solution of DDEs

I have the following problem:
I solved DDEs with some parameters which have two solutions: r and s.
Now, I need to calculate and plot (in Manipulate) the following integral:

And there are problems. Which function to use to calculate the integral? NIntegrate doesn’t work because “x in limit is unexpected”. If I try to use Integrate I obtain the following error:

Integral and error estimates are 0 on all integration subregions. Try
increasing the value of the MinRecursion option. If value of integral
may be 0, specify a finite value for the AccuracyGoal option.

But Integrate does not have MinRecursion and AccuracyGoal options. Also, I’m not sure if I put functions r and s in the integral correctly (see code below, these functions in the final Manipulate should be of h as in prob function?). How do I make Manipulate work?

CODE:
Auxiliary functions:

hill[x_, n_] := x^n/(1 + x^n)
input = 0.4
steadypar[n_, hillweight_] := 
 Sort[{r, s} /. 
   NSolve[{input - r + s * hillweight * hill[r * s, n] == 0, 
     input - s + r * hillweight * hill[r * s, n] == 0, r >= 0, 
     s >= 0}, {r, s}, Reals]]

My DDEs solver:

modelsolver[σ_, τ_, hillweight_, end_, sigend_] := 
 NDSolve[{r'[t] == 
    input  + Piecewise[{{σ, t < sigend}, {0, t >= 0.5}}] - 
     r [t - τ] + s[t] * hillweight * hill[r[t] * s[t], 2], 
   s'[t] == 
    input - s[t - τ] + r[t] * hillweight * hill[r[t] * s[t], 2], 
   r[t /; t <= 0] == First[steadypar[2, hillweight]][[1]], 
   s[t /; t <= 0] == First[steadypar[2, hillweight]][[2]]}, {r, 
   s}, {t, 0, end}]

Function that should calculate integral:

probtest[x_, β_, r_, s_] := 1/(1 + E^(-β * Integrate[r - s, {h, 0, x}]))

Alternative version:

prob[x_?NumericQ, β_, r_, s_] := 
 1/(1 + E^(-β * 
       NIntegrate[r - s, {h, 0, x}, MinRecursion -> 100, 
        AccuracyGoal -> 100]))

Final Manipulate:

Simulation = 
 Manipulate[
  Plot[Evaluate[({probtest[x, β, r[h], s[h]], 
       1 - probtest[x, β, r[h], s[h]]} /. 
      modelsolver[σ, τ, hillweight, end, 0.5][[1]])], {x, 
    0, end}, PlotRange -> {0, 1}, MaxRecursion -> 15 ], {τ, 0, 
   2, 0.1}, {σ, 0 , 4, 0.01}, {hillweight, 0.85, 1, 
   0.01}, {end, 10, 200, 15}, {β, 1, 1.3, 0.05}]

modelsolver[σ_, τ_, hillweight_, end_, sigend_] := NDSolve[
    {
    r'[t] == input  + Piecewise[{{σ, t < sigend}, {0, t >= 0.5}}] - r[t - τ] + s[t] * hillweight * hill[r[t] * s[t], 2], 
    s'[t] == input - s[t - τ] + r[t] * hillweight * hill[r[t] * s[t], 2], 
    r[t /; t <= 0] == First[steadypar[2, hillweight]][[1]], 
    s[t /; t <= 0] == First[steadypar[2, hillweight]][[2]],
    int'[t] == r[t] - s[t],
    int[t /; t<=0] == 0
    },
    {r, s, int},
    {t, 0, end}
]

sol = modelsolver[2, 1, .85, 10, .5];
Plot[int[t] /. sol, {t, 0, 10}, PlotRange->All]