def METHOD_NAME(egf_doses, vem_doses, readout='absolute'):
    erk = np.zeros((len(egf_doses), len(vem_doses)))
    ras = np.zeros((len(egf_doses), len(vem_doses)))
    for i, ed in enumerate(egf_doses):
        yobs0, y0 = simulate_untreated(model, ts, ed)
        for j, vd in enumerate(vem_doses):
            y0_last = y0[-1]
            y0_last[model.observables['Vem_free'].species[0]] = vd
            yobs_ss, yobs_tc = get_steady_state(model, y0_last)
            if readout == 'absolute':
                erk[i,j] = yobs_ss['ERK_p']
                ras[i,j] = yobs_ss['RAS_active']
            elif readout == 'foldchange':
                erk[i,j] = yobs_ss['ERK_p'] / yobs0['ERK_p'][-1]
                ras[i,j] = yobs_ss['RAS_active'] / yobs0['RAS_active'][-1]
            elif readout == 'ss_min_diff':
                erk[i,j] = yobs_ss['ERK_p'] - yobs_tc['ERK_p'].min()
                ras[i,j] = yobs_ss['RAS_active']
    return erk, ras