def _pypsa_load_timeseries(network, timesteps, mode=None):
    """
    Time series in PyPSA compatible format for load instances

    Parameters
    ----------
    network : Network
        The eDisGo grid topology model overall container
    timesteps : array_like
        Timesteps is an array-like object with entries of type
        :pandas:`pandas.Timestamp<timestamp>` specifying which time steps
        to export to pypsa representation and use in power flow analysis.
    mode : str, optional
        Specifically retrieve load time series for MV or LV grid level or both.
        Either choose 'mv' or 'lv'.
        Defaults to None, which returns both timeseries for MV and LV in a
        single DataFrame.

    Returns
    -------
    :pandas:`pandas.DataFrame<dataframe>`
        Time series table in PyPSA format
    """
    mv_load_timeseries_p = []
    mv_load_timeseries_q = []
    lv_load_timeseries_p = []
    lv_load_timeseries_q = []

    # add MV grid loads
    if mode is 'mv' or mode is None:
        for load in network.mv_grid.graph.nodes_by_attribute('load'):
            mv_load_timeseries_q.append(load.pypsa_timeseries('q').rename(
                repr(load)).to_frame().loc[timesteps])
            mv_load_timeseries_p.append(load.pypsa_timeseries('p').rename(
                repr(load)).to_frame().loc[timesteps])
        if mode is 'mv':
            lv_load_timeseries_p, lv_load_timeseries_q = \
                _pypsa_load_timeseries_aggregated_at_lv_station(
                    network, timesteps)

    # add LV grid's loads
    if mode is 'lv' or mode is None:
        for lv_grid in network.mv_grid.lv_grids:
            for load in lv_grid.graph.nodes_by_attribute('load'):
                lv_load_timeseries_q.append(load.pypsa_timeseries('q').rename(
                    repr(load)).to_frame().loc[timesteps])
                lv_load_timeseries_p.append(load.pypsa_timeseries('p').rename(
                    repr(load)).to_frame().loc[timesteps])

    load_df_p = pd.concat(mv_load_timeseries_p + lv_load_timeseries_p, axis=1)
    load_df_q = pd.concat(mv_load_timeseries_q + lv_load_timeseries_q, axis=1)

    return load_df_p, load_df_q