So you can use ParametricNDSolve, but it's not very well suited to being nested since it doesn't return just a simple interpolation. If you were still interested in knowing how you can convert an Integral into a Parametric NDSolve statement, please consider this integral as an example:
Integrate[Sin[a x], {x, 0, 1}]
This can be easily solved symbolically, returning a function of a, but let's just assume it can't. If that were the case, we could transform it into this ParametricNDSolveValue statement:
ParametricNDSolveValue[{output'[x] == Sin[a x], output[0] == 0}, output[1], {x, 0, 1}, {a}]
This gives the same curve, just numerically.
Anyway, what you probably really want to do is build an Interpolation on each iteration. There's smarter ways of doing this, but I'll try a rough example with equidistant, nonadaptive sampling.
Consider this simple integral equation which I'll use as an example because it is really simple. I don't know that much about integral equations or I would come up with something more interesting that was roughly as simple
F[a] == Integral[F[x a], {x, 0, 1}]
We want to make a function, I'll call it iterate, which takes a function and then runs:
NIntegrate[f[x a], {x, 0, 1}]
For many values of "a" and then makes an interpolation out of those values. Use Table to create a list of values for different values of a. It can be defined like this:
iterate[f_] := Interpolation@
Table[{a, NIntegrate[f[x a], {x, 0, 1}]}, {a, 0, 1, 0.01}]
"iterate" is a function that returns a function defined on the interval 0 to 1
We can now run iterate over and over. I'll just choose Sin as the seed function:
it = iterate[Sin]
Now run this a couple of times
it = iterate[it]; Plot[it[t], {t, 0, 1}]
You can see that it approaches a function equal to 0 everywhere, which is boring, but a solution to the intergral equation above.
