How to model species interactions (interaction matrix) for inverse problem modelling with SciML software

Question

I have a probably simple question that I haven't found an answer to. I'm building ODEs for a competitive Lotka Volterra model in order to perform inverse modelling. I first built a 3 species dummy model spelling out each parameter one by one which works really well but I am having trouble expanding to a greater number of species using broadcasting. I suspect that it is because I'm attempting to reference parameters from a matrix and an array since my attempt at broadcasting a logistic growth model with only an array of growth rates was successful.

As far as I can tell I've managed to correctly specify the model. I can't figure out how to feed in starting parameters to my problem

#LotkaVolterra competition model (self-interacting terms and K set to 1)

using ModelingToolkit, LinearAlgebra

n = 3 #number of species considered
u0 = repeat([0.01], n)
tsteps = [0:1:100;]
tspan = (0, last(tsteps))

###Species interaction matrix (we ignore self interaction)
m=ones(n,n)
m[diagind(m)] .= 0.0


@variables t, N(..)[1:n]
@parameters begin
    ρ[1:n,], [bounds=(0, Inf), tunable=true]
    α[1:n, 1:n], [bounds=(0, 1), tunable=true]
end

D = Differential(t)

vec_eqs = D.(N(t)) .~ (1 .-N(t) .-m .*α *N(t)).* ρ .*N(t)

@named LV_mod = ODESystem(vec_eqs)


#LV_prob = ODEProblem(structural_simplify(LV_mod), u0, tspan, ???)

I assume that there should be some way to correctly integrate values from the initial guesses for the ρ array and α matrix (here with values of 0.25 everywhere).

ρ_ini= repeat([0.25], n)
α_ini =repeat([0.25], n, n)
α_ini[diagind(α_ini)] .= 0.0

Any help greatly appreciated! ...also happy for any feedback on code optimization, I am very much a beginner to SciML and Julia.

Update:

Since I was not considering the diagonal of my species interaction matrix I just coded the diagonal as my species-dependant growth rates. I can now input my parameters in the problem as a single matrix and then do inverse modelling. After optimization the model fits the data well. My updated parameters are then spit out as a vector. Can I just use the reshape function on the result.u vector to recover my initial matrix?

I would still love to figure out how to best specify the parameters (or models?) to accommodate for more complex/flexible parameter structures...and recover updated parameters in their original format if possible.

Sorry for the rambling and likely misuse of terms!