Is it possible to solve the following system of integrodifferential equations in Mathematica?

Consider a system of equations
$$
begin{cases}F(E,t) = frac{partial f(E,t)}{partial t}-H(T_{gamma})Efrac{partial f(E,t)}{partial E} -I(f(E,t),E,T_{gamma})=0, \ G(t) = frac{partial T_{gamma}}{partial t} + H(T_{gamma})T_{gamma}+int limits_{0}^{infty} dE I(f(E,t),E,T_{gamma})=0,end{cases}
$$

where
$$
I(f(E,t),E,T_{gamma}) equivint limits_{Gamma(E)} dE’ G_{1}(E’,E,f(E,t),f(E’,t),T_{gamma})
$$

is some integral and $H$ is a known function and $Gamma(E)$ is the domain of integration.

Could you please tell me whether it is possible to solve this system in Mathematica? As far as I know, NDSolve does not work with integrodifferential systems of equations, but maybe some scheme exists?

These equations in Mathematica code, Equationfn(En, t), EquationTg(t) are given below (there is some piece of code that defines functions):

ME = 0.5;
etaBPlanck = 6.1*10^-10;
ng(T_) = 2*Zeta(3)/Pi^2 T^3;
rhoeTemp(T_) := 
 4/(2*Pi^2)
   NIntegrate((Ee^2 Sqrt(Ee^2 - ME^2))/(
   Exp(Ee/T) + 1), {Ee, ME, Infinity})
peTemp(T_) := 
 4/(6*Pi^2)
   NIntegrate((Ee^2 - ME^2)^(3/2)/(Exp(Ee/T) + 1), {Ee, ME, Infinity})
rhoe(T_) = 
 Quiet(10^Interpolation(
     Join(Table({10^T, Log10(etaBPlanck*ng(10^T)*ME)}, {T, 
        Log10(10^-8), Log10(0.0049), 0.1}), 
      Table({10^T, 
        Log10(etaBPlanck*ng(10^T)*ME + rhoeTemp(10^T))}, {T, 
        Log10(0.005), Log10(19.99), 0.01}), 
      Table({T, Log10(7/8*4*Pi^2/30*T^4)}, {T, 20, 300, 0.5})), 
     InterpolationOrder -> 1)(T))
pe(T_) = Quiet(
  10^Interpolation(
     Join(Table({10^T, -90}, {T, Log10(10^-8), Log10(0.0049), 0.1}), 
      Table({10^T, Log10(10^-99 + peTemp(10^T))}, {T, Log10(0.005), 
        Log10(19.99), 0.01}), 
      Table({T, Log10(1/3 7/8*4*Pi^2/30*T^4)}, {T, 20, 300, 0.5})), 
     InterpolationOrder -> 1)(T))
gStare(T_) := rhoe(T)/(Pi^2/30 T^4)
DrhoeDT(T_) = D(rhoe(T), T);
(*Neutrino and photon*)
rhon(T_) = 7/8*6*Pi^2/30*T^4;
rhona(T_) = 7/8*2*Pi^2/30 T^4;
rhog(T_) = 2*Pi^2/30 T^4;
DrhogDT(T_) = D(rhog(T), T);

GFvalue = 1.16*10^-11;
glValue = 1/2 + 0.22;
grValue = -(1/2) + 0.22;
hbar = 6.582119*10^-25;
sToGeVminusOne = 1/hbar;
sToMeVminusOne = sToGeVminusOne*10^-3;
CollisionIntegralnnToee(En_, 
  t_) := -fn(En, t)*
   NIntegrate(
    sToMeVminusOne (
     8 En E2^3 GFvalue^2 (glValue^2 + grValue^2) fn(E2, t))/(
     9 (Pi)^3), {E2, 0, Infinity}) + (
  16 E^(-(En/Tg(t))) En *sToMeVminusOne*
   GFvalue^2 (glValue^2 + grValue^2) Tg(t)^4)/(3 (Pi)^3)
CollisionIntegralneTone(En_, 
  t_) := -((
   64 En*sToMeVminusOne*
    GFvalue^2 (glValue^2 + grValue^2) Tg(t)^4 fn(En, t))/(
   3 (Pi)^3)) + 
  NIntegrate( 
   sToMeVminusOne/(En^2 (Pi)^3) 8 E^(-(En/Tg(t)))
     GFvalue^2 (glValue^2 + grValue^2) Tg(
     t)^2 (-En^2 (E3^2 + 2 E3 Tg(t) - 
         2 (-1 + E^(E3/Tg(t))) Tg(t)^2) - 
      2 En Tg(t) (E3^2 + 2 E3 Tg(t) - 
         2 (-1 + E^(E3/Tg(t))) Tg(t)^2) + 
      2 Tg(t)^2 ((-1 + E^(E3/Tg(t))) E3^2 - 
         2 (1 + E^(E3/Tg(t))) E3 Tg(t) + 
         4 (-1 + E^(E3/Tg(t))) Tg(t)^2)) fn(E3, t), {E3, 0, En}) + 
  NIntegrate(
   sToMeVminusOne/(En^2 (Pi)^3) 8 E^(-(En/Tg(t)))
     GFvalue^2 (glValue^2 + grValue^2) Tg(
     t)^2 (2 (-1 + E^(En/Tg(t))) Tg(
        t)^2 (E3^2 + 2 E3 Tg(t) + 4 Tg(t)^2) - 
      En^2 (E3^2 + 2 E3 Tg(t) - 2 (-1 + E^(En/Tg(t))) Tg(t)^2) - 
      2 En Tg(t) (E3^2 + 2 E3 Tg(t) + 
         2 (1 + E^(En/Tg(t))) Tg(t)^2)) fn(E3, t), {E3, En, Infinity}) 
H(t_) = sToMeVminusOne (1.66 Sqrt(10.75) Tg(t)^2)/(1.2*10^22);
Equationfn(En_, t_) := ((D(fn(En, tt), tt)) /. {tt -> t}) - 
  H(Tg(t)) En*((D(fn(en, t), en)) /. {en -> En}) - 
  CollisionIntegralnnToee(En, t) - CollisionIntegralneTone(En, t)==0
EquationTg(t_) := ((D(Tg(tt), tt)) /. {tt -> t}) + 
  1/(DrhoeDT(Tg(t)) + 
    DrhogDT(Tg(t))) (4*H(Tg(t))*rhog(Tg(t)) + 
     3*H(Tg(t))*(rhoe(Tg(t)) + pe(Tg(t))) + 
     NIntegrate(
      CollisionIntegralnnToee(En, t) En^3/(2*Pi^2), {En, 0, 
       Infinity}) + 
     NIntegrate(
      CollisionIntegralneTone(En, t) En^3/(2*Pi^2), {En, 0, 
       Infinity}))==0