Lotka Volterra estimate parameters from experimental data

[quantumCircuit]

Homework Statement

I am trying model Gause's experiment about the paramecium for interspecific competition according to Lotka-Volterra model, but I am stuck at estimating the parameters of the system of differential equations.

Relevant Equations

The system of equations is the following:
$\dot{x} = x(L-Ax-By)$
$\dot{y} = y(M-Cx-Dy)$

Namely, in the system, I have obtained the value of parameters L, M, A and D, because I treat the other organism as equal to zero, i.e., it doesn't exist, but I am struggling about the values of B and C, that are coupled with the product of x and y. Can anyone help me how to obtain those values. I am programming in Python, and I used curve_fit to fit the logistic equation.

 

[BvU]

Hello @quantumCircuit , !

to fit the logistic equation.

So I suppose you have set(s) of data of x(t) and y(t) ?
Maybe even with different $x(0)$ and $y(0)$ ?

And from your exposé I understand you did a preliminary fit to
$$\begin{align*} \dot{x} &= x(L-Ax-By) \\ \dot{y} &= y(M-Cx-Dy) \end{align*} $$with $B = C = 0$ ? And that gives an unsatisfactory fit ?

 

[pasmith]

Fitting $x$ and $y$ to the logistic equation separately is not going to work; the non-linear interaction makes that impossible, since $x$ and $y$ will likely exhibit behaviour which two uncoupled logistic systems cannot produce (such as periodic orbits or convergence to a fixed point which is not on the coordinate axes).

For the 1D problem you can solve the ODE analytically and then use your favorite regression algorithm to estimate the paramters; but for the 2D problem I don't think it's possible to find $(x(t),y(t))$ other than by solving the ODE numerically from given initial conditions using estimated parameter values, comparing the result to your data, and working out how to change the paramters to improve the fit.

You can place bounds on the parameter values by looking at the qualitative behaviour of your data: does it tend to a solution where both species are extinct, where only one is extinct, to a fixed point where neither species is extinct, or does it converge to a periodic solution?

 

[quantumCircuit]

Fitting $x$ and $y$ to the logistic equation separately is not going to work; the non-linear interaction makes that impossible, since $x$ and $y$ will likely exhibit behaviour which two uncoupled logistic systems cannot produce (such as periodic orbits or convergence to a fixed point which is not on the coordinate axes).

For the 1D problem you can solve the ODE analytically and then use your favorite regression algorithm to estimate the paramters; but for the 2D problem I don't think it's possible to find $(x(t),y(t))$ other than by solving the ODE numerically from given initial conditions using estimated parameter values, comparing the result to your data, and working out how to change the paramters to improve the fit.

You can place bounds on the parameter values by looking at the qualitative behaviour of your data: does it tend to a solution where both species are extinct, where only one is extinct, to a fixed point where neither species is extinct, or does it converge to a periodic solution?

Yes, this is exactly what I am trying to do, but the only problem being is that I don't have estimated parameters, so I am trying to do the estimation. I have already tried to estimate them with setting x=0 and y=0, as I said before, and also I have tried to simultaneously use least square method on the equations, but I didn't obtain satisfactory results, not even remotely. Thus, my question remains the same: How does one estimate parameters of already given experimental data, and known system of differential equations, specifically Lotka-Volterra model for competition?

 

[BvU]

Can you perhaps be bothered to answer the questions in post #2 ?
Do you understand where they come from ?
Do you understand what is mentioned by @pasmith in post #3 ?

Can you solve the equations ?

(I mean the other ones, sorry)
I would start with an ODE solver and minimize $\sum \text{residuals}^2$,

But first of all I would look at the data; time plots & phase-spasce plots.

 

[quantumCircuit]

Can you perhaps be bothered to answer the questions in post #2 ?
Do you understand where they come from ?
Do you understand what is mentioned by @pasmith in post #3 ?

Can you solve the equations ?

(I mean the other ones, sorry)
I would start with an ODE solver and minimize $\sum \text{residuals}^2$,

But first of all I would look at the data; time plots & phase-spasce plots.

About post #2, no, it doesn't give satisfactory result because the system is coupled, and to answer the question in #3, yes, on of the species is going extinct, and the other thrives. I don't know how to solve a coupled differential system, as far as I have searched on the Internet, I couldn't find exact, analytic solution.

Could you suggest a starting point for a solver to minimize the residuals for a system of nonlinear differential equations?

 

[BvU]

I googled 'python solving ode'. Plenty hits, e.g. (indirect):

https://docs.scipy.org/doc/scipy/re...rate.solve_ivp.html#scipy.integrate.solve_ivp

The scenario: The solver gives you $\vec y(t)$ the calculated arrays $x$ and $y$.
You calculate the residuals and wrap that calculation in a function that returns $\sum \text{residuals}^2$.
Then let optimize.minimize mimize that function

I've done that many times, but with fortran (python beginner). There should be examples in python, maybe even completely worked out. Just google a bit.

Perhaps this junky site ?

 

[quantumCircuit]

I googled 'python solving ode'. Plenty hits, e.g. (indirect):

https://docs.scipy.org/doc/scipy/re...rate.solve_ivp.html#scipy.integrate.solve_ivp

The scenario: The solver gives you $\vec y(t)$ the calculated arrays $x$ and $y$.
You calculate the residuals and wrap that calculation in a function that returns $\sum \text{residuals}^2$.
Then let optimize.minimize mimize that function

I've done that many times, but with fortran (python beginner). There should be examples in python, maybe even completely worked out. Just google a bit.

Perhaps this junky site ?

Thank you for the reply, but I would like to ask you for a little clarification here, as this is something I am completely unfamiliar with. I cannot get the minimize function working. So far I have done this. I have the experimental data for x(t), and y(t), and I have values for the parameters. I managed to find somewhat satisfactory values from trial and error starting from an educated guess. Still, I am not satisfied with those results.

What would be my objective function to use the residuals? Do I need two objective functions, one for x(t), and the other for y(t)? I calculated x(t) (experimental) - x(t) (estimated), squared these differences, and then summed them. But, I don't really know how to proceed from here on. How to use the minimize function from Python with this objective function?

 

[BvU]

Do I need two objective functions, one for x(t), and the other for y(t)?

You need one function, $$\sum_i \Bigl (x(t_i)_\text{exp}-x(t_i)_\text{calc} \Bigr)^2 + \sum_i \Bigl (y(t_i)_\text{exp}-y(t_i)_\text{calc} \Bigr)^2$$

Re python: did you look at a few examples ? They call it 'fun' there
Do you understand that the argument for the function is an array with L M A B C D E ?
Initial guesses go in $x0$, optimized values should come back in res.x

I have the experimental data for x(t), and y(t), and I have values for the parameters.

Can you post something, like a plot of $\ x(t_i)_\text{exp},\ x(t_i)_\text{calc}, \ y(t_i)_\text{exp}, \ y(t_i)_\text{calc}\$ for the initial guesses ?

Do you have only one experiment, or many ? All with $x_0, y_0 = 0$ ?

 

[BvU]

I managed to write the objective function, and it is working, and I am getting specific values for the parameters, (1) but, it appears that my calculated values somehow don't agree with the minimize function. I calculated the data by using the Euler direct method (2) with the estimated values for the parameters, and when I send that to the minimize function (3) , I get complete nonsense. I suppose you would like to see some code, or plots, but I just want to ask you right now what could possibly go wrong? (4)

I prefer the public thread: ohers might help and/or correct us when we go astray !

With my limited Python knowledge it's hard to guess what goes wrong and where. On top of that my telepathic capabilities are non-existent .
Using libraries like scipy is all nice and well, but getting control is cumbersome. Best advice is to follow the calculations step by step and check every time.

Can you post something, like a plot of $\ x(t_i)_\text{exp},\ x(t_i)_\text{calc}, \ y(t_i)_\text{exp}, \ y(t_i)_\text{calc}\$ for the initial guesses ?

And the corresponding objective function value ?

(1) The objective function should have the parameters LMABCD as input, not output ?
Ideally you should let it also simultaneously plot $\ x(t_i)_\text{exp},\ x(t_i)_\text{calc}, \ y(t_i)_\text{exp}, \ y(t_i)_\text{calc}\$ and print the return value so you can manually vary LMABCD and check the behaviour.

Do your measurements have errors, are there outliers, etc ?

(2) Is that an option in scipy.integrate.solve_ivp or did you write your own integrator ?

(3) Sounds unlogical: you are at a minimum and let the optimizer start from there ?

(4) That's easy: according to the well-established Murphy laws, anything that can go wrong will . Bugs don't print messages like "hey this is wrong"; instead they have to be rooted out one by one.

 

[quantumCircuit]

I prefer the public thread: ohers might help and/or correct us when we go astray !

With my limited Python knowledge it's hard to guess what goes wrong and where. On top of that my telepathic capabilities are non-existent .
Using libraries like scipy is all nice and well, but getting control is cumbersome. Best advice is to follow the calculations step by step and check every time.

And the corresponding objective function value ?

(1) The objective function should have the parameters LMABCD as input, not output ?
Ideally you should let it also simultaneously plot $\ x(t_i)_\text{exp},\ x(t_i)_\text{calc}, \ y(t_i)_\text{exp}, \ y(t_i)_\text{calc}\$ and print the return value so you can manually vary LMABCD and check the behaviour.

Do your measurements have errors, are there outliers, etc ?

(2) Is that an option in scipy.integrate.solve_ivp or did you write your own integrator ?

(3) Sounds unlogical: you are at a minimum and let the optimizer start from there ?

(4) That's easy: according to the well-established Murphy laws, anything that can go wrong will . Bugs don't print messages like "hey this is wrong"; instead they have to be rooted out one by one.

I think the problem is that the function minimize from SciPy doesn't take the data I am sending it. It just tries to find the minimum of a 4-value function with six parameters, but it doesn't calculate the residuals of the experimental and calculated data at all.

 

[BvU]

but it doesn't calculate

You are supposed to do that calculation yourself in a function fun that is callable by minimize.

Take an example and replace fun by your own written function.

Can you post your code, so I or someone else can help improve ?
( post using [code=python]blbla[/code] )

 

[quantumCircuit]

You are supposed to do that calculation yourself in a function fun that is callable by minimize.

Take an example and replace fun by your own written function.

Can you post your code, so I or someone else can help improve ?
( post using [code=python]blbla[/code] )

Here is the code I used so far:

def objective_function(z): 
    return (np.sum((z[0] - z[1])**2) + np.sum((z[2] - z[3])**2))/2

z = np.array([x_exp, x_calc, y_exp, y_calc])

opt = minimize(objective_function, x0=params)
opt.x

 

[BvU]

I don't see objective_function re-calculate $x_\text{calc}$ and $y_\text{calc}$ with the updated values for the parameters ?

It seems to me array z doesn't change between calls to objective_function

 

[quantumCircuit]

I don't see objective_function re-calculate $x_\text{calc}$ and $y_\text{calc}$ with the updated values for the parameters ?

It seems to me array z doesn't change between calls to objective_function

Yes, but shouldn't that be a job of the minimize function?

 

[BvU]

No. Minimize receives the sum of squares and modifies the parameters LMABCD until a minimum is achieved.

Look at the examples.