How to plot a vector field using Plots.jl and function with DynamicalSystems.jl signature

I am trying to convert the python code for a tutorial on Numerical bifurcation diagrams and bistable systems to Julia; however, I am struggling to make a vector field/quiver/streamplot using Julia. How can I do this? A full example of my attempt is included below.

using DynamicalSystems
using Plots
using UnPack
using LaTeXStrings

"""
NEED HELP HERE
"""
function plot_flow!(myplot, α, β, xspace = 0:0.1:2, yspace = 0:0.1:2)
    flow = [
         cellular_switch_rule(
            [ui, vi],
            Dict(:α => α, :β => β), 
            nothing) 
        for ui in xspace, vi in yspace]
    # not sure what to do about the indexing here?
    us = [flow[i, j][1] for i in 1:length(xspace), j in 1:length(yspace)]
    vs = [flow[i, j][2] for i in 1:length(xspace), j in 1:length(yspace)]
    quiver!(myplot, xspace, yspace, quiver=(us, vs))
end 

"""
    cellular_switch_rule(y, p, t)

Biochemical rate equation formulation of gene expression for E-coli. 
Compliant with DynamicalSystems.jl expected function signature.
"""
function cellular_switch_rule(y, p, t)
    @unpack α, β = p
    u, v = y
    udot = α/(1 + v^β) - u
    vdot = α/(1 + u^β) - v
    return SVector(udot, vdot)
end 

cellular_switch(u0, p) = CoupledODEs(cellular_switch_rule, u0, p)

u_nullcline(v, α, β) = α/(1 + v^β)
v_nullcline(u, α, β) = α/(1 + u^β)

function plot_isoclines!(myplot, α, β, uspace = 0:0.1:2, vspace = 0:0.1:2)
    n_u = [u_nullcline(vi, α, β) for vi in vspace]
    n_v = [v_nullcline(ui, α, β) for ui in uspace]
    plot!(myplot, uspace, n_v, 
        linstyle = :dash,
        label = L"\frac{dv}{dt} = 0 \Leftrightarrow n_v(u)")
    plot!(myplot, n_u, vspace, 
        linestyle = :dash,
        label = L"\frac{du}{dt} = 0 \Leftrightarrow n_u(v)")
end 

# Attempt to plot nullclines + vector field
α = 1
β = 2

myplot = plot(
    xlabel = "u", ylabel = "v", 
    xlim = (0, 2), ylim = (0, 2),
    legend = :right)

plot_isoclines!(myplot, α, β)
plot_flow!(myplot, α, β)

Unfortunately, the quiver function from Plots.jl while functional upon correction of the way I handled indexing and through the creation of a proper mesh Julia style, see below,

# julia style meshgrid
xxs = [x for x in xspace, y in yspace]; 
yys = [y for x in xspace, y in space]

the resulting Plots.quiver plot does not scale the arrows very well, and it turns out to be far too messy to read. Since I had trouble using the PyPlot.streamplot function, I decied to use CairoMakie.streamplot, which is really intuitive to use and the working code is shown below:

using DynamicalSystems: CoupledODEs, trajectory, SVector
using CairoMakie; cm = CairoMakie
using UnPack
using LaTeXStrings

# parameters
α = 1
β = 2

# global parameter dict for use with streamplot
P = Dict(:α => α, :β => β)

# Functions needed for streamplot
function cellular_switch_point2f(y, P) 
    Point2f(cellular_switch_rule(y, P, nothing))
end 
    
flow_field(x) = cellular_switch_point2f(x, P)

"""
    cellular_switch_rule(u, p, t)

Biochemical rate equation formulation of gene expression for E-coli. 

# References

[1] : http://www.normalesup.org/~doulcier/teaching/modeling/bistable_systems.html

[2] : Gardner, T., Cantor, C. & Collins, J. Construction of a genetic 
toggle switch in Escherichia coli. Nature 403, 339–342 (2000). 
https://doi.org/10.1038/35002131
"""
function cellular_switch_rule(y, p, t) 
    @unpack α, β = p
    u, v = y
    udot = α/(1 + v^β) - u
    vdot = α/(1 + u^β) - v
    return SVector(udot, vdot)
end 
    
cellular_switch(u0, p) = CoupledODEs(cellular_switch_rule, u0, p)

u_nullcline(v, α, β) = α/(1 + v^β)
v_nullcline(u, α, β) = α/(1 + u^β)
    
function plot_isoclines!(ax, α, β, uspace = 0:0.1:2, vspace = 0:0.1:2)
    n_u = [u_nullcline(vi, α, β) for vi in vspace]
    n_v = [v_nullcline(ui, α, β) for ui in uspace]
    cm.lines!(
        ax, 
        uspace, 
        n_v,
        linestyle = :dash,
        linewidth = 3,
        label = L"\frac{dv}{dt} = 0 \Leftrightarrow n_v(u)")
    cm.lines!(
        ax, 
        n_u, 
        vspace,
        linestyle = :dash,
        linewidth = 3, 
        label = L"\frac{du}{dt} = 0 \Leftrightarrow n_u(v)")
end 
    
fig = cm.Figure()
ax = Axis(
    fig[1, 1],
    xlabel = "u",
    ylabel = "v",
    title = "Nullclines & Flow Field of E-coli Model"
)
cm.xlims!(ax, (0, 2))
cm.ylims!(ax, (0, 2))
        
plot_isoclines!(ax, α, β)
cm.streamplot!(ax, flow_field, 0..2, 0..2, colormap = :grays, alpha=0.2)
cm.axislegend(ax)
fig