The auto bifurcation diagram seems doesn't catch all multi-steady states?

[ZiyiiiLiu]

Hi, I used the auto to plot the bifurcation diagram of my ode set, as the figures attached that there is no multi-stability when S1 =1.5 based on the bifurcation diagram, which should mean that varying the initial conditions the steady states should keep the same.
But when I changed the initial condition, the level of each gene is different.
For example, S1=1.5, when initial condition :ics_1 = {'Y': 0.01, 'A': 1.1, 'N': 1.5, 'Z': 0.1, 'C': 1, 'V': 0.1, 'H': 2}, the level of A is 0.6 which fits the bifurcation diagram, but when ics_1 = {'Y': 0.01, 'A': 3.1, 'N': 1.5, 'Z': 0.1, 'C': 1, 'V': 0.1, 'H': 2}, the level of A is 3.7, these results shows there should be at least 2 states when S1=1.5 which didn't fit the bifurcation diagram?

I attached my code here:
model_sclc_2 = """
k_0Y = 0.1; s0_Y = 0.5; K_Y_S1 = 1; n_Y_S1 = 4; K_AY = 4; n_AY = 5;
k_0A = 0.4; K_A = 4.5; K_AA = 1.7; n_AA = 2.7; K_NA = 1.5; n_NA = 3.3;
K_HA = 1.5; n_HA = 3.3; r_dA = 1.3; k_0N = 1.09; K_N = 2.5; K_NN = 0.2;
n_NN = 5.84; K_AN = 3.77; n_AN = 5.61; r_dN = 0.24; K_A_S1 = 3.4; n_A_S1 = 3.85;
s0_N = 0.46; K_N_S1 = 1.94; n_N_S1 = 4.52; S1 = 1.2;
k_0Z = 0.001; s0_Z = 0.1; K_Z_N = 1; n_Z_N = 3; K_Z_S = 1; n_Z_S = 1; r_dZ = 2;
k_0C = 1; s0_C = 1; K_S1_C = 1; n_S1_C = 1; r_dC = 1;
k_0V = 1; s0_V = 1; K_V_S1 = 1; n_V_S1 = 1; K_AV = 1; n_AV = 1; r_dV = 1;
k_0H = 1; s0_H = 1; K_AH = 1; n_AH = 1; K_NH = 1; n_NH = 1; r_dH = 1;

-> Y; k_0Y + s0_Y * (S1 / K_Y_S1)^n_Y_S1 / (1 + (S1 / K_Y_S1)^n_Y_S1) * (1 / (1 + (A / K_AY)^n_AY)) - 0.2 * Y
-> A; k_0A + K_A * (A / K_AA)^n_AA / (1 + (A / K_AA)^n_AA) * ((1 / (1 + (N / K_NA)^n_NA)) + (1 / (1 + (S1 / K_A_S1)^n_A_S1)) * (1 / (1 + (H / K_HA)^n_HA))) - r_dA * A
-> N; k_0N + (K_N * (N / K_NN)^n_NN / (1 + (N / K_NN)^n_NN) + (S1 + s0_N) / (1 + (S1 / K_N_S1)^n_N_S1)) * (1 / (1 + (A / K_AN)^n_AN) + 1 / (1 + (S1 / K_A_S1)^n_A_S1)) - r_dN * N
-> Z; k_0Z + s0_Z * (N / K_Z_N)^n_Z_N / (1 + (S1 / K_Z_N)^n_Z_N) * (1 / (1 + (S1 / K_Z_S)^n_Z_S)) - r_dZ * Z
-> C; k_0C + s0_C / (1 + (S1 / K_S1_C)^n_S1_C) - r_dC * C
-> V; k_0V + s0_V * (S1 / (1 + (S1 / K_V_S1)^n_V_S1)) * (1 / (1 + (A / K_AV)^n_AV)) - r_dV * V
-> H; k_0H + s0_H * (A / K_AH)^n_AH / (1 + (A / K_AH)^n_AH) * (1 / (1 + (N / K_NH)^n_NH)) - r_dH * H

"""
r = te.loada(model_sclc_2)

ics_1 = {'Y': 0.01, 'A':1.1, 'N': 1.5, 'Z': 0.1, 'C': 1, 'V': 0.1, 'H': 2}
for var, val in ics_1.items():
setattr(r, var, val)

pars = {
'k_0Y': 0.1, 's0_Y': 0.8, 'K_Y_S1': 1.6, 'n_Y_S1': 5, 'K_AY': 5, 'n_AY': 5,
'k_0A': 0.826, 'K_A': 3.82, 'K_AA': 1.8, 'n_AA': 3, 'K_NA': 3, 'n_NA': 5,
'K_HA': 1, 'n_HA': 3, 'r_dA': 1.59, 'k_0N': 0.01, 'K_N': 0.7, 'K_NN': 2.5,
'n_NN': 5, 'K_AN': 3.5, 'n_AN': 5, 'r_dN': 0.99, 'K_A_S1': 1.5, 'n_A_S1': 5,
's0_N': 0.9, 'K_N_S1': 1.5, 'n_N_S1': 5, 'S1': 1.5, 'k_0Z': 0.0249,
's0_Z': 0.8237, 'K_Z_N': 3.5, 'n_Z_N': 4, 'K_Z_S': 1.4568, 'n_Z_S': 4.6874,
'r_dZ': 2, 'k_0C': 0.1665, 's0_C': 0.5962, 'K_S1_C': 0.33, 'n_S1_C': 4,
'r_dC': 0.55, 'k_0V': 0.1, 's0_V': 1, 'K_V_S1': 3, 'n_V_S1': 4, 'K_AV': 1,
'n_AV': 4, 'r_dV': 1, 'k_0H': 0.1, 's0_H': 1, 'K_AH': 3, 'n_AH': 3,
'K_NH': 1, 'n_NH': 4, 'r_dH': 1
}
for k in pars:
setattr(r, k, pars[k])
m = r.simulate(0, 1000, 5000)

time = m[:, 0] # Time data
Y = m[:, 1] # 'Y' variable data
A = m[:, 2] # 'A' variable data
N = m[:, 3] # 'N' variable data
Z = m[:, 4] # 'Z' variable data
C = m[:, 5] # 'C' variable data
V = m[:, 6] # 'V' variable data
H = m[:, 7] # 'H' variable data

plt.figure(figsize=(16,3))

plt.plot(time, A, label='A')
plt.plot(time, N, label='N ')
plt.plot(time, Z, label='Z ')
plt.plot(time,H, label='H ')
plt.plot(time,Y, label='Y ')
ics_title = ', '.join([f'{key}={val}' for key, val in ics_1.items()])
plt.xlabel('Time')
plt.ylabel('Concentration')
#plt.title('Simulation Results of SCLC Model')
plt.title(f'Simulation Results of SCLC Model with Initial Conditions: {ics_title}')
plt.legend(loc='best')
plt.grid(True)
plt.show()

auto = Plugin("tel_auto2000")
auto.setProperty("SBML", r.getCurrentSBML())
auto.setProperty("ScanDirection", "Positive")
auto.setProperty("PrincipalContinuationParameter", "S1")
auto.setProperty("PCPLowerBound", 0)
auto.setProperty("PCPUpperBound", 5)
auto.setProperty("NMX", 10000)
auto.execute()
print("Summary: ", auto.BifurcationSummary)
pts = auto.BifurcationPoints
lbls = auto.BifurcationLabels
biData = auto.BifurcationData
biData.plotBifurcationDiagram(pts, lbls)