Is there a more efficient way to implement the introduced model in GurobiPy?

I want to implement the following model in GurobiPy which is a minimal example of the problem I am confronted with:

As the data I work with has a large number of set elements for IJK, JKL, and KLM, I need an efficient implementation for summing over my x variables that satisfy the index mapping ((j,k) in IJK and (j,k,l) in JKL and (k,l,m) in KLM). Iterating over all elements of sets IJK, JKL, and KLM seems to become inefficient as the cardinality of these sets increases.

The best performing implementation I came up with is as follows:

import gurobipy as gpy
def model(I, IJK, JKL, KLM):
    model = gpy.Model()
    x_list = [(i, j, k, l, m)
              for (i, j, k) in IJK.select('*', '*', '*')
              for (j, k, l) in JKL.select(j, k, '*')
              for (k, l, m) in KLM.select(k, l, '*')]

    x = model.addVars(x_list, name='x')

    model.setObjective(1, gpy.GRB.MINIMIZE)

    model.addConstrs(
        (gpy.quicksum(x[k] for k in x_list if k[0] == i) >= 0 for i in I),
        'ei'
    )

    model.update()

Here is code to produce random data for the model:

import numpy as np
import pandas as pd
import gurobipy as gpy

def create_random_data(n, m):
    I = [f'i{x}' for x in range(1, n + 1)]
    J = [f'j{x}' for x in range(1, m + 1)]
    K = [f'k{x}' for x in range(1, m + 1)]
    L = [f'l{x}' for x in range(1, m + 1)]
    M = [f'm{x}' for x in range(1, m + 1)]

    ijk = pd.DataFrame(np.random.binomial(1, 0.05, size=(len(I)*len(J)*len(K))),
                       index=pd.MultiIndex.from_product(
                           [I, J, K], names=['i', 'j', 'k']),
                       columns=['value']).reset_index()
    jkl = pd.DataFrame(np.random.binomial(1, 0.05, size=(len(J)*len(K)*len(L))),
                       index=pd.MultiIndex.from_product(
                           [J, K, L], names=['j', 'k', 'l']),
                       columns=['value']).reset_index()
    klm = pd.DataFrame(np.random.binomial(1, 0.05, size=(len(K)*len(L)*len(M))),
                       index=pd.MultiIndex.from_product(
                           [K, L, M], names=['k', 'l', 'm']),
                       columns=['value']).reset_index()

    IJK = gpy.tuplelist([tuple(x) for x in ijk.loc[ijk['value'] == 1][['i', 'j', 'k']].to_dict('split')['data']])
    JKL = gpy.tuplelist([tuple(x) for x in jkl.loc[jkl['value'] == 1][['j', 'k', 'l']].to_dict('split')['data']])
    KLM = gpy.tuplelist([tuple(x) for x in klm.loc[klm['value'] == 1][['k', 'l', 'm']].to_dict('split')['data']])
    
    return I, J, K, L, M, IJK, JKL, KLM

Is there a more efficient way to implement the introduced model in GurobiPy?

I, J, K, L, M, IJK, JKL, KLM = create_random_data(300, 20)
%timeit -r 7 -n 10 model(I, IJK, JKL, KLM)

297 ms ± 66.9 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)