How to use scipy.integrate.odeint(func, y0, t, args=(),...) when there are some constraints among the parameters of the function "func"?

Question

Here is my attempt to solve a system of ordinary differential equations given some initial conditions and constraints. While I have figured out how to solve the ODE using odeint from Scipy package without the constraints, I do not know how to incorporate the existence of two constraints on the parameters of the function (which are commented out in the code). I would truly appreciate your help.

Note: any undefined variable in the code is actually some fixed scalar.

---

import numpy as np
from scipy.integrate import odeint

#initial conditions
y0, u0 = ([2.1E19, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [2.1E19, 0, 0, 0])


def lhs(y):
    """This function extracts all 13 components of the solution y"""
    return y[0], y[1], y[2], y[3], y[4], y[5], y[6], y[7], y[8], y[9], y[10], y[11], y[12]


def eqns(time_val, y, f_func, g_func, h_func, p_func, diff_eq_solutions_num):
    "This is the system of equations that needs to be solved numerically"
    var1, var2, var3, var4, var5, var6, var7, var8, var9, var10, var11, var12, var13 = lhs(y)
    #var1
    r0 = -(K2 * var11 * var1) - ((K4+K5) * var11 * var1) + (0.5 * (K6+K7) * var2**2) + (
        2*K12*var11*var8) - (K11 * var1 * (var5+var6+var7)) + (Qa*var2*var9) + (Qb*var3*var1) + (
            Qc*var4*var1) + (A4*var7) + (A5*var6) - (
                (K14+K15) * var11 * var1) - (f_func(time_val) * var1) - (
                    g_func(time_val) * var1) - (G1+G2+G3+G4) * h_func(time_val)
    #var2
    r1 = (A2*var3) - ((K6+K7) * var2**2) - (K11 * var2 * (var5+var6+var7)) + (
        K3 * f1 * var11 *
        (var5+var6+var7)) + (K14*var11*var1) - (Qa*var9*var2) + G2 * f1 * h_func(time_val)
    #var3
    r2 = (A1*var4) - (A2*var3) + (0.5 * K6 * var2**2) - (K11 * var3 * (var5+var6+var7)) - (
        Qb*var3*var1) + (K10 * var1 * var9**2) + (K15*var11*var1) + (
            K3 * f2 * var11 * (var5+var6+var7)) + G2 * f2 * h_func(time_val)
    #var4
    r3 = (K3 * f3 * var11 * (var5+var6+var7)) - (A1*var4) - (Qc*var4*var1) + (K7*var2*var2) - (
        K11 * var4 * (var5+var6+var7)) + (f_func(time_val) * var1) + G2 * f3 * h_func(time_val)
    #var5
    r4 = (K2*fi1*var11*var1) - (K3*var11*var5) - (K8*var11*var5) - (K11*var1*var5) + (A5*var6) + (
        A4*var7) + (Qi1*var7*var1) + (Qi2*var7*NO2) + (fi1 * g_func(time_val) *
                                                       var1) + G1 * fi1 * h_func(time_val)
    #var6
    r5 = (K2*fi2*var11*var1) - (K3*var11*var6) - (A5*var6) - (K8*var11*var6) - (K11*var1*var5) + (
        fi2 * g_func(time_val) * var1) + G1 * fi3 * h_func(time_val)
    #var7
    r6 = (K2*fi3*var11*var1) - (K3*var11*var7) - (A4*var7) - (Qi1*var7*var1) - (K8*var11*var7) - (
        Qi2*var7*NO2) + (fi3 * g_func(time_val) * var1) + G1 * fi2 * h_func(time_val)
    #var8
    r7 = (K11 * (var5+var6+var7) * (var1+var2+var3+var4)) - (K12*var11*var8)
    #var9
    r8 = (K4*var11*var1) + (2*K5*var11*var1) + (2 * K8 * var11 * (var5+var6+var7)) - (
        K9*var11*var9) - (K10 * var1 * var9**2) + (K13*var11*var10) + G4 * h_func(time_val)
    #var10
    r9 = (K4*var11*var1) + (K9*var11*var9) - (K13*var11*var10) + G3 * h_func(time_val)
    r10 = -(K8 * var11 * (var5+var6+var7)) + (K9*var11*var9) - (K13*var11*var10) - (
        K12*var11*var8) - (K3 * var11 *
                           (var5+var6+var7)) + p_func(time_val) + G5 * h_func(time_val) + (
                               K4*var11*var1) + (K2*var11*var1)
    r11 = (A1*var4) - (4 * np.pi * var12 / 3E-11)
    r12 = (A4*var7) - (4 * np.pi * var13 / 3E-11)
    #These are the constraints I would like to incorporate into system of differential equations
    # 2.1E19 = var1 + var2 + var3 + var4 + var5 + var6 + var7 + var8 + var9 + var10
    # 0 = var5 + var6 + var7 + var8 + var10 - var11
    return [r0, r1, r2, r3, r4, r5, r6, r7, r8, r9, r10, r11, r12]


def integrator(soln, f_func, g_func, p_func, h_func, i, j, diff_eq_solutions_num):
    args = (f_func[i][j], g_func[i][j], p_func[i][j], h_func[i][j], diff_eq_solutions_num)
    t = np.linspace(0, 2E-6, 100000)
    soln[i][j] = odeint(eqns, y0, t, args=args, tfirst=True, printmessg=True, ixpr=True)