diff --git a/tvb_library/tvb/simulator/models/stefanescu_jirsa.py b/tvb_library/tvb/simulator/models/stefanescu_jirsa.py
index cac8e956ac..45eff9d60d 100644
--- a/tvb_library/tvb/simulator/models/stefanescu_jirsa.py
+++ b/tvb_library/tvb/simulator/models/stefanescu_jirsa.py
@@ -33,6 +33,7 @@
 from scipy.stats import norm as scipy_stats_norm
 from .base import Model
 from tvb.basic.neotraits.api import NArray, Final, List, Range
+from numba import njit
 
 
 class ReducedSetBase(Model):
@@ -193,10 +194,7 @@ class ReducedSetFitzHughNagumo(ReducedSetBase):
 
     def dfun(self, state_variables, coupling, local_coupling=0.0):
         r"""
-
-
         The system's equations for the i-th mode at node q are:
-
         .. math::
                 \dot{\xi}_{i}    &=  c\left(\xi_i-e_i\frac{\xi_{i}^3}{3} -\eta_{i}\right)
                                   + K_{11}\left[\sum_{k=1}^{o} A_{ik}\xi_k-\xi_i\right]
@@ -211,9 +209,7 @@ def dfun(self, state_variables, coupling, local_coupling=0.0):
                                   + \left[\sum_{k=1}^{o} W_{\zeta}\cdot\xi_{kr}\right] \\
                                  & \\
                 \dot{\beta}_i    &= \frac{1}{c}\left(\alpha_i-b\beta_i+n_i\right)
-
         """
-
         xi = state_variables[0, :]
         eta = state_variables[1, :]
         alpha = state_variables[2, :]
@@ -221,7 +217,6 @@ def dfun(self, state_variables, coupling, local_coupling=0.0):
         derivative = numpy.empty_like(state_variables)
         # sum the activity from the modes
         c_0 = coupling[0, :].sum(axis=1)[:, numpy.newaxis]
-
         # TODO: generalize coupling variables to a matrix form
         # c_1 = coupling[1, :] # this cv represents alpha
 
@@ -240,15 +235,80 @@ def dfun(self, state_variables, coupling, local_coupling=0.0):
 
         return derivative
 
+    @njit
+    def dfun_numba(self, state_variables, coupling, local_coupling=0.0):
+        r"""
+
+
+        The system's equations for the i-th mode at node q are:
+
+        .. math::
+                \dot{\xi}_{i}    &=  c\left(\xi_i-e_i\frac{\xi_{i}^3}{3} -\eta_{i}\right)
+                                  + K_{11}\left[\sum_{k=1}^{o} A_{ik}\xi_k-\xi_i\right]
+                                  - K_{12}\left[\sum_{k =1}^{o} B_{i k}\alpha_k-\xi_i\right] + cIE_i                       \\
+                                 &\, + \left[\sum_{k=1}^{o} \mathbf{\Gamma}(\xi_{kq}, \xi_{kr}, u_{qr})\right]
+                                  +  \left[\sum_{k=1}^{o} W_{\zeta}\cdot\xi_{kr} \right] \\
+                \dot{\eta}_i     &= \frac{1}{c}\left(\xi_i-b\eta_i+m_i\right)                                              \\
+                & \\
+                \dot{\alpha}_i   &= c\left(\alpha_i-f_i\frac{\alpha_i^3}{3}-\beta_i\right)
+                                  + K_{21}\left[\sum_{k=1}^{o} C_{ik}\xi_i-\alpha_i\right] + cII_i                          \\
+                                 & \, + \left[\sum_{k=1}^{o} \mathbf{\Gamma}(\xi_{kq}, \xi_{kr}, u_{qr})\right]
+                                  + \left[\sum_{k=1}^{o} W_{\zeta}\cdot\xi_{kr}\right] \\
+                                 & \\
+                \dot{\beta}_i    &= \frac{1}{c}\left(\alpha_i-b\beta_i+n_i\right)
+
+        """
+
+        xi = state_variables[0, :]
+        eta = state_variables[1, :]
+        alpha = state_variables[2, :]
+        beta = state_variables[3, :]
+        derivative = numpy.empty_like(state_variables)
+        # sum the activity from the modes
+        c_0 = coupling[0, :].sum(axis=1)[:, numpy.newaxis]
+
+        # TODO: generalize coupling variables to a matrix form
+        # c_1 = coupling[1, :] # this cv represents alpha
+
+        N = len(xi)
+        
+        # Compute numpy.dot(xi, self.Aik) manually
+        xi_Aik = numpy.zeros(N)
+        for i in range(N):
+            xi_Aik[i] = sum(xi[j] * self.Aik[j, i] for j in range(N))
+
+        # Compute numpy.dot(alpha, self.Bik) manually
+        alpha_Bik = numpy.zeros(N)
+        for i in range(N):
+            alpha_Bik[i] = sum(alpha[j] * self.Bik[j, i] for j in range(N))
+
+        # Compute numpy.dot(xi, self.Cik) manually
+        xi_Cik = numpy.zeros(N)
+        for i in range(N):
+            xi_Cik[i] = sum(xi[j] * self.Cik[j, i] for j in range(N))
+
+        for i in range(N):
+            derivative[0][i] = (self.tau * (xi[i] - self.e_i * xi[i] ** 3 / 3.0 - eta[i]) +
+                                self.K11 * (xi_Aik[i] - xi[i]) -
+                                self.K12 * (alpha_Bik[i] - xi[i]) +
+                                self.tau * (self.IE_i + c_0 + local_coupling * xi[i]))
+
+            derivative[1][i] = (xi[i] - self.b * eta[i] + self.m_i) / self.tau
+
+            derivative[2][i] = (self.tau * (alpha[i] - self.f_i * alpha[i] ** 3 / 3.0 - beta[i]) +
+                                self.K21 * (xi_Cik[i] - alpha[i]) +
+                                self.tau * (self.II_i + c_0 + local_coupling * xi[i]))
+
+            derivative[3][i] = (alpha[i] - self.b * beta[i] + self.n_i) / self.tau
+            return derivative
+
     def update_derived_parameters(self):
         """
         Calculate coefficients for the Reduced FitzHugh-Nagumo oscillator based
         neural field model. Specifically, this method implements equations for
         calculating coefficients found in the supplemental material of
         [SJ_2008]_.
-
         Include equations here...
-
         """
         newaxis = numpy.newaxis
         trapz = scipy_integrate_trapz
@@ -281,7 +341,6 @@ def update_derived_parameters(self):
         g2 = norm.pdf(Zu)
         G1 = numpy.tile(g1, (self.number_of_modes, 1))
         G2 = numpy.tile(g2, (self.number_of_modes, 1))
-
         cV = numpy.conj(V)
         cU = numpy.conj(U)
 
@@ -302,6 +361,115 @@ def update_derived_parameters(self):
 
         self.m_i = (self.a * intcVdZ).T
         self.n_i = (self.a * intcUdZ).T
+
+    @njit
+    def trapz_integrate_numba(self,y, x):
+        """Compute the trapezoidal integral explicitly."""
+        result = 0.0
+        for i in range(len(x) - 1):
+            result += 0.5 * (x[i+1] - x[i]) * (y[i+1] + y[i])
+        return result
+
+    @njit
+    def update_derived_parameters_numba(self):
+        """
+        Calculate coefficients for the Reduced FitzHugh-Nagumo oscillator based
+        neural field model. Specifically, this method implements equations for
+        calculating coefficients found in the supplemental material of
+        [SJ_2008]_.
+
+        Include equations here...
+
+        """
+
+        stepu = 1.0 / (self.nu + 2 - 1)
+        stepv = 1.0 / (self.nv + 2 - 1)
+
+        norm = scipy_stats_norm(loc=self.mu, scale=self.sigma)
+        
+        # Generate Zu and Zv explicitly
+        Zu = numpy.zeros(self.nu)
+        Zv = numpy.zeros(self.nv)
+        
+        for i in range(self.nu):
+            Zu[i] = norm.ppf((i + 1) * stepu)
+        
+        for i in range(self.nv):
+            Zv[i] = norm.ppf((i + 1) * stepv)
+        
+        # Initialize U and V matrices
+        V = numpy.zeros((self.number_of_modes, self.nv))
+        U = numpy.zeros((self.number_of_modes, self.nu))
+
+        nv_per_mode = self.nv // self.number_of_modes
+        nu_per_mode = self.nu // self.number_of_modes
+
+        # Assign ones in blocks
+        for i in range(self.number_of_modes):
+            for j in range(i * nv_per_mode, (i + 1) * nv_per_mode):
+                V[i, j] = 1.0
+            for j in range(i * nu_per_mode, (i + 1) * nu_per_mode):
+                U[i, j] = 1.0
+
+        # Normalize the modes using explicit integration
+        for i in range(self.number_of_modes):
+            V[i, :] /= numpy.sqrt(self.trapz_integrate(V[i, :] * V[i, :], Zv))
+            U[i, :] /= numpy.sqrt(self.trapz_integrate(U[i, :] * U[i, :], Zu))
+
+        # Compute normal PDFs
+        g1 = numpy.zeros(self.nv)
+        g2 = numpy.zeros(self.nu)
+
+        for i in range(self.nv):
+            g1[i] = norm.pdf(Zv[i])
+        
+        for i in range(self.nu):
+            g2[i] = norm.pdf(Zu[i])
+
+        # Preallocate matrices
+        self.Aik = numpy.zeros((self.number_of_modes, self.number_of_modes))
+        self.Bik = numpy.zeros((self.number_of_modes, self.number_of_modes))
+        self.Cik = numpy.zeros((self.number_of_modes, self.number_of_modes))
+
+        self.e_i = numpy.zeros((1, self.number_of_modes))
+        self.f_i = numpy.zeros((1, self.number_of_modes))
+        self.IE_i = numpy.zeros((1, self.number_of_modes))
+        self.II_i = numpy.zeros((1, self.number_of_modes))
+        self.m_i = numpy.zeros((self.number_of_modes, 1))
+        self.n_i = numpy.zeros((self.number_of_modes, 1))
+
+        # Compute conjugates
+        cV = numpy.conj(V)
+        cU = numpy.conj(U)
+
+        # Compute integrations
+        intcVdZ = numpy.zeros((self.number_of_modes, 1))
+        intG1VdZ = numpy.zeros((1, self.number_of_modes))
+        intcUdZ = numpy.zeros((self.number_of_modes, 1))
+
+        for i in range(self.number_of_modes):
+            intcVdZ[i, 0] = self.trapz_integrate(cV[i, :], Zv)
+            intG1VdZ[0, i] = self.trapz_integrate(g1 * V[i, :], Zv)
+            intcUdZ[i, 0] = self.trapz_integrate(cU[i, :], Zu)
+
+        # Compute Aik, Bik, Cik using explicit loops
+        for i in range(self.number_of_modes):
+            for j in range(self.number_of_modes):
+                self.Aik[i, j] = intcVdZ[i, 0] * intG1VdZ[0, j]
+                self.Bik[i, j] = intcVdZ[i, 0] * self.trapz_integrate(g2 * U[j, :], Zu)
+                self.Cik[i, j] = intcUdZ[i, 0] * intG1VdZ[0, j]
+
+        # Compute e_i, f_i, IE_i, II_i
+        for i in range(self.number_of_modes):
+            self.e_i[0, i] = self.trapz_integrate(cV[i, :] * (V[i, :] ** 3), Zv)
+            self.f_i[0, i] = self.trapz_integrate(cU[i, :] * (U[i, :] ** 3), Zu)
+            self.IE_i[0, i] = self.trapz_integrate(Zv * cV[i, :], Zv)
+            self.II_i[0, i] = self.trapz_integrate(Zu * cU[i, :], Zu)
+
+        # Compute m_i and n_i
+        for i in range(self.number_of_modes):
+            self.m_i[i, 0] = self.a * intcVdZ[i, 0]
+            self.n_i[i, 0] = self.a * intcUdZ[i, 0]
         # import pdb; pdb.set_trace()
 
 
@@ -488,7 +656,6 @@ class ReducedSetHindmarshRose(ReducedSetBase):
     def dfun(self, state_variables, coupling, local_coupling=0.0):
         r"""
         The equations of the population model for i-th mode at node q are:
-
         .. math::
                 \dot{\xi}_i     &=  \eta_i-a_i\xi_i^3 + b_i\xi_i^2- \tau_i
                                  + K_{11} \left[\sum_{k=1}^{o} A_{ik} \xi_k - \xi_i \right]
@@ -507,9 +674,7 @@ def dfun(self, state_variables, coupling, local_coupling=0.0):
                 & \\
                 \dot{\beta}_i   &= h_i - p_i \alpha_i^2 - \beta_i \\
                 \dot{\gamma}_i  &= rs \alpha_i - r \gamma_i - n_i
-
         """
-
         xi = state_variables[0, :]
         eta = state_variables[1, :]
         tau = state_variables[2, :]
@@ -546,9 +711,7 @@ def update_derived_parameters(self, corrected_d_p=True):
         of Hindmarsh-Rose oscillators. Specifically, this method implements
         equations for calculating coefficients found in the supplemental
         material of [SJ_2008]_.
-
         Include equations here...
-
         """
 
         newaxis = numpy.newaxis
@@ -565,7 +728,6 @@ def update_derived_parameters(self, corrected_d_p=True):
         # Define the modes
         V = numpy.zeros((self.number_of_modes, self.nv))
         U = numpy.zeros((self.number_of_modes, self.nu))
-
         nv_per_mode = self.nv // self.number_of_modes
         nu_per_mode = self.nu // self.number_of_modes
 
@@ -613,6 +775,185 @@ def update_derived_parameters(self, corrected_d_p=True):
         else:
             # typo in the original paper by RS & VJ, kept for comparison purposes.
             self.d_i = (self.d * intcVdI).T
+
+    @njit
+    def dfun_numba(self, state_variables, coupling, local_coupling=0.0):
+        r"""
+        The equations of the population model for i-th mode at node q are:
+
+        .. math::
+                \dot{\xi}_i     &=  \eta_i-a_i\xi_i^3 + b_i\xi_i^2- \tau_i
+                                 + K_{11} \left[\sum_{k=1}^{o} A_{ik} \xi_k - \xi_i \right]
+                                 - K_{12} \left[\sum_{k=1}^{o} B_{ik} \alpha_k - \xi_i\right] + IE_i \\
+                                &\, + \left[\sum_{k=1}^{o} \mathbf{\Gamma}(\xi_{kq}, \xi_{kr}, u_{qr})\right]
+                                 + \left[\sum_{k=1}^{o} W_{\zeta}\cdot\xi_{kr} \right] \\
+                & \\
+                \dot{\eta}_i    &=  c_i-d_i\xi_i^2 -\tau_i \\
+                & \\
+                \dot{\tau}_i    &=  rs\xi_i - r\tau_i -m_i \\
+                & \\
+                \dot{\alpha}_i  &=  \beta_i - e_i \alpha_i^3 + f_i \alpha_i^2 - \gamma_i
+                                 + K_{21} \left[\sum_{k=1}^{o} C_{ik} \xi_k - \alpha_i \right] + II_i \\
+                                &\, +\left[\sum_{k=1}^{o}\mathbf{\Gamma}(\xi_{kq}, \xi_{kr}, u_{qr})\right]
+                                 + \left[\sum_{k=1}^{o}W_{\zeta}\cdot\xi_{kr}\right] \\
+                & \\
+                \dot{\beta}_i   &= h_i - p_i \alpha_i^2 - \beta_i \\
+                \dot{\gamma}_i  &= rs \alpha_i - r \gamma_i - n_i
+
+        """
+
+        xi = state_variables[0, :]
+        eta = state_variables[1, :]
+        tau = state_variables[2, :]
+        alpha = state_variables[3, :]
+        beta = state_variables[4, :]
+        gamma = state_variables[5, :]
+        derivative = numpy.empty_like(state_variables)
+        
+        c_0 = numpy.sum(coupling[0, :], axis=1)[:, numpy.newaxis]
+
+        for i in range(len(xi)):
+            sum_Aik_xi = 0.0
+            sum_Bik_alpha = 0.0
+            sum_Cik_xi = 0.0
+            
+            for k in range(len(xi)):
+                sum_Aik_xi += xi[k] * self.A_ik[i, k]
+                sum_Bik_alpha += alpha[k] * self.B_ik[i, k]
+                sum_Cik_xi += xi[k] * self.C_ik[i, k]
+            
+            derivative[0, i] = (eta[i] - self.a_i[i] * xi[i] ** 3 + self.b_i[i] * xi[i] ** 2 - tau[i] +
+                                self.K11 * (sum_Aik_xi - xi[i]) -
+                                self.K12 * (sum_Bik_alpha - xi[i]) + self.IE_i[i] +
+                                c_0[i, 0] + local_coupling * xi[i])
+            
+            derivative[1, i] = self.c_i[i] - self.d_i[i] * xi[i] ** 2 - eta[i]
+            
+            derivative[2, i] = self.r * self.s * xi[i] - self.r * tau[i] - self.m_i[i]
+            
+            derivative[3, i] = (beta[i] - self.e_i[i] * alpha[i] ** 3 + self.f_i[i] * alpha[i] ** 2 - gamma[i] +
+                                self.K21 * (sum_Cik_xi - alpha[i]) + self.II_i[i] +
+                                c_0[i, 0] + local_coupling * xi[i])
+            
+            derivative[4, i] = self.h_i[i] - self.p_i[i] * alpha[i] ** 2 - beta[i]
+            
+            derivative[5, i] = self.r * self.s * alpha[i] - self.r * gamma[i] - self.n_i[i]
+        
+        return derivative
+
+    @njit
+    def trapezoidal_integral(self,y, x):
+        """Manually compute the trapezoidal integral."""
+        integral = 0.0
+        for i in range(len(x) - 1):
+            integral += 0.5 * (x[i + 1] - x[i]) * (y[i + 1] + y[i])
+        return integral
+
+    @njit
+    def normalize_modes(self,modes, I_values):
+        """Normalize the modes manually."""
+        num_modes, size = modes.shape
+        for i in range(num_modes):
+            norm_factor = numpy.sqrt(self.trapezoidal_integral(modes[i] * modes[i], I_values))
+            if norm_factor > 0:
+                modes[i] /= norm_factor
+        return modes
+
+    @njit
+    def manual_dot(self,vec1, mat):
+        """Manually computes the dot product of a vector with a matrix."""
+        result = numpy.zeros(mat.shape[1])
+        for j in range(mat.shape[1]):
+            for i in range(mat.shape[0]):
+                result[j] += vec1[i] * mat[i, j]
+        return result
+
+    @njit
+    def update_derived_parameters_numba(self, corrected_d_p=True):
+        """
+        Calculate coefficients for the neural field model based on a Reduced set
+        of Hindmarsh-Rose oscillators. Specifically, this method implements
+        equations for calculating coefficients found in the supplemental
+        material of [SJ_2008]_.
+
+        Include equations here...
+
+        """
+
+        stepu = 1.0 / (self.nu + 2 - 1)
+        stepv = 1.0 / (self.nv + 2 - 1)
+
+        norm = self.norm_dist
+        Iu = numpy.zeros(self.nu)
+        Iv = numpy.zeros(self.nv)
+        for i in range(self.nu):
+            Iu[i] = norm.ppf((i + 1) * stepu)
+        for i in range(self.nv):
+            Iv[i] = norm.ppf((i + 1) * stepv)
+
+        # Define the modes
+        V = numpy.zeros((self.number_of_modes, self.nv))
+        U = numpy.zeros((self.number_of_modes, self.nu))
+
+        nv_per_mode = self.nv // self.number_of_modes
+        nu_per_mode = self.nu // self.number_of_modes
+
+        for i in range(self.number_of_modes):
+            for j in range(i * nv_per_mode, (i + 1) * nv_per_mode):
+                if j < self.nv:
+                    V[i, j] = 1.0
+            for j in range(i * nu_per_mode, (i + 1) * nu_per_mode):
+                if j < self.nu:
+                    U[i, j] = 1.0
+
+        # Normalise the modes
+        V = self.normalize_modes(V, Iv)
+        U = self.normalize_modes(U, Iu)
+
+        # Compute Gaussian PDFs
+        g1 = numpy.zeros_like(Iv)
+        g2 = numpy.zeros_like(Iu)
+        for i in range(len(Iv)):
+            g1[i] = norm.pdf(Iv[i])
+        for i in range(len(Iu)):
+            g2[i] = norm.pdf(Iu[i])
+
+        # Compute conjugates
+        cV = numpy.conj(V)
+        cU = numpy.conj(U)
+
+        # Compute integrals manually
+        intcVdI = numpy.zeros((self.number_of_modes, 1))
+        intcUdI = numpy.zeros((self.number_of_modes, 1))
+        intG1VdI = numpy.zeros((1, self.number_of_modes))
+
+        for i in range(self.number_of_modes):
+            intcVdI[i, 0] = self.trapezoidal_integral(cV[i], Iv)
+            intcUdI[i, 0] = self.trapezoidal_integral(cU[i], Iu)
+            intG1VdI[0, i] = self.trapezoidal_integral(g1 * V[i], Iv)
+
+        # Compute coefficients
+        self.A_ik = self.manual_dot(intcVdI, intG1VdI).T
+        self.B_ik = self.manual_dot(intcVdI, self.trapezoidal_integral(g2 * U, Iu)[numpy.newaxis, :])
+        self.C_ik = self.manual_dot(intcUdI, intG1VdI).T
+
+        self.a_i = self.a * self.trapezoidal_integral(cV * V ** 3, Iv)[numpy.newaxis, :]
+        self.e_i = self.a * self.trapezoidal_integral(cU * U ** 3, Iu)[numpy.newaxis, :]
+        self.b_i = self.b * self.trapezoidal_integral(cV * V ** 2, Iv)[numpy.newaxis, :]
+        self.f_i = self.b * self.trapezoidal_integral(cU * U ** 2, Iu)[numpy.newaxis, :]
+        self.c_i = (self.c * intcVdI).T
+        self.h_i = (self.c * intcUdI).T
+
+        self.IE_i = self.trapezoidal_integral(Iv * cV, Iv)[numpy.newaxis, :]
+        self.II_i = self.trapezoidal_integral(Iu * cU, Iu)[numpy.newaxis, :]
+
+        if corrected_d_p:
+            # correction identified by Shrey Dutta & Arpan Bannerjee, confirmed by RS
+            self.d_i = self.d * self.trapezoidal_integral(cV * V ** 2, Iv)[numpy.newaxis, :]
+            self.p_i = self.d * self.trapezoidal_integral(cU * U ** 2, Iu)[numpy.newaxis, :]
+        else:
+            # typo in the original paper by RS & VJ, kept for comparison purposes.
+            self.d_i = (self.d * intcVdI).T
             self.p_i = (self.d * intcUdI).T
 
         self.m_i = (self.r * self.s * self.xo * intcVdI).T