David,

Sorry that I missed the attachment.

This whole issue is exactly why I tried to do this in Mathematica to see if there were differences and there are. I wanted the most accurate answer I could get and Mathematica gives me that.

That Euler notebook in Mathematica looks complicated. I have some studying to do there!

I am satisfied my original question has been answered thanks to both Bill and you!

Scott

David,

This is beautiful work and what I have been after for about 5 days now! To me it seems clear from your two simulations using Euler's method that one never knows when the time increment is small enough to ensure accuracy. One also runs into computational time problems in Excel with really, really small time increments. I guess I would have to go to 10 million rows in Excel to get the accuracy I need. Thank you again for all the time you put into solving this conundrum for me.

Do you mind sending me your notebook where you implemented Euler's method? I would like to know how to do that method in Mathematica for comparison purposes.

Scott

Hi Scott,

I implemented Euler's method in Mathematica using a function called NestList. It repeatedly applies a function to the last result, starting with some initial condition. I used your equations for the derivatives and your time step. The original result that both Bill and I have using NDSolve is this:

The result using Euler's method is this:

You can see that they are quite similar. And they converge on the same concentrations. I suspect the difference arises from the fixed time steps being to large in the early part of the simulation, resulting in under damping. Mathematica has adaptive time stepping to prevent that.

I don't know why the Excel Euler's method produces different results, but I would trust NDSolve over the perhaps-too-simple Euler's method.

Attachments:

By Euler.nb

David,

To save you time looking for the corrected truncated Excel sheet, here it is again.

Scott

Attachments:

Oscillating equa...xlsx

I'll have a look. Although I find reading Excel is a bit difficult. I am looking at doing Eulers method with Mathematica using your time step.

With that fine a time step and reasonably well behaved equations I doubt that the problem is in the accuracy of Euler's method.

David and Bill,

Here is a Youtube link showing the method (Euler's) I used to solve those three differential equations in Excel. https://www.youtube.com/watch?v=vXWb1Kspx8s I think you surmised this already, Bill, when you deconstructed my Excel equations. Euler's method is crude and slow, but I always thought if the dt was small enough, it approximated derivatives pretty accurately (hence my small time increment of 0.000025 min. and my 1 million row Excel sheet).

Scott

Thank you for all your expertise and help! I am glad it is fun for you. Not so much for me yet. Hopefully, someday I will get there. I will take a look at the notebook in an hour or two. I am really glad you engaged in this help request.

Scott

Bill,

David is correct. All the K parameters should be multiplied by their respective protein concentrations in all the differential equations.

Scott

I am guessing I still don't have this correct because the plots do not match. If you can find any mistakes that I have made then please point them out. I'm trying to reproduce your very first example.

If you can think of any little experiment we can try on simpler made-up systems then maybe that would confirm that simpler examples do match between Excel and Mathematica.

If I look at the equations for each row of your Excel sheet then I think I see

t+dt,((BX*MX/(1+KZ*zm)-xm/TX)*dt+xm)/dt,((BY*MY/(1+KX*xm)-ym/TY)*dt+ym)/dt,((BY*MY/(1+KX*ym)-zm/TY)*dt+zm)/dt

I think that is what I am trying to reproduce in Mathematica, with the translation from discrete to differential equation. But it still seems like I am misunderstanding something.

Thank you for the text. That was invaluable.

In textbook Figure 8-83

d[X]/dt=betax mx 1/(1+(Ky[Y])^hy-[X]/taux
d[Y]/dt=betay my 1/(1+(Kx[X])^hx-[Y]/tauy

Mathematica

x'[t]==betax mx / (1+ky[y[t]]^hy)-x[t]/taux
y'[t]==betay my / (1+kx[x[t]]^hy)-y[t]/tauy

Excel Sheet image of 3 equations in upper right

dx/dt=betax mx/(1+Kz[z])-[x]/taux
dy/dt=betay my/(1+Kx[x])-[y]/tauy
dz/dt=betaz mz/(1+Ky[y])-[z]/tauz

Mathematica

x'[t]==betax mx/(1+kz[z[t]])-x[t]/taux
y'[t]==betay my/(1+kx[x[t]])-y[t]/tauy
z'[t]==betaz mz/(1+ky[y[t]])-z[t]/tauz

I am guessing in your Excel sheet that you chose

betax=1
betay=1
betaz=1
mx=1
my=1
mz=1
taux=2
tauy=1.5
tauz=2

Now I need to understand your K functions. Is

Kx[z]=10^8*z?

That doesn't seem right. Likewise for Ky and Kz

It looks like he is describing a variety of different situations that can arise and I'm not certain yet just what of those you are trying to do in your model

If I can understand those last three and I haven't made any other mistakes then I think I'm close to pasting this into Mathematica and getting a plot.

David,

Could you also show me with one of my parameters (pick the one you want) how to add a slider bar with the Manipulate command? I believe you used pretty fancy coding of my equations, and I tried for about 90 min. to make the Manipulate function work and could not. I was able to get the 12 parameter sliders added, but somehow they are not linked to the calculation and do nothing. When I trained with the tutorials in Mathematica, the Manipulate function was always presented as Manipulate[Plot[Evaluate.....] and your notebook is written pretty differently from that example so I cannot figure out how to get Manipulate to work. If you give me one example, I can add the other 11 parameters since i want to be able to adjust the 12 parameters that go into this calculation to see what effect each has on the oscillation.

Sorry for the bother again!

Scott

I've tried and can't get my hands on a copy of your textbook at the moment.

If you could copy the differential equations from the text into a message here and include the values of the parameters that you are using at the moment then I will try to translate that system into Mathematica and plot the solutions and see how close my plots are to your posted Excel plots. I might be able to generate your entire Excel sheet from the example you gave, but I'm not sure I can do that correctly.

Hopefully we are going to get you what you need. Thanks for your patience and assistance with this.

Oops -- you were typing this while I was responding above. ;-)

David,

Due to my inexperience with Mathematica, the last notebook I sent you did not evaluate correctly. In the newly attached notebook, i do see all proteins overlapping with one another, so that is pretty convincing evidence to me that Mathematica is working correctly. i wish i was a better mathematician because I still do not understand why Excel is not giving the correct curves. The only think I can think of is that 0.000025 is still not a small enough time increment to allow Excel to simulate the differential equations accurately.

In any event, I think I will believe the Mathematica curves more than the curves in Excel. Thank you again!

Scott

Attachments:

677785769-DK all...nb

David,

Thank you for the input. It makes sense what you are saying. I went ahead and made all the parameters the same for each type of parameter for each of the three proteins. Under those conditions, all the proteins should reach the same steady-state level. You will see in the attached notebook they do not, so I am still puzzled by what is going on in Mathematica. When I simulated in Excel, all three proteins reached the same-steady level, which should be the case since all proteins are being made and broken down at the same rate. I wanted to do this comparison to try to figure out why there is disagreement between Excel and Mathematica.

Any further feedback would be appreciated.

Scott

Attachments:

677785769-DK SK ...nb

Hi Scott,

While x and z do have the same lifetimes, the coupling of equations breaks the symmetry. x'[t] depends on z[t], but z'[t] depends an y[t]. The equations with the parameters substituted look like this:

If you set the derivatives equal to zero, and solve for the three concentrations, they are not the same. (That is what DFixedPoint does. Its new in 14.1.) Maybe I have made an error in the equations or parameters.

Regarding time stepping, I did not specify time step parameters to Mathematica. It has very clever algorithms for doing that.

Hi Scott,

I've also looked at this. When I use your equations and parameter values from your spreadsheet graphic, I see this: What is most puzzling is that your graph appears to show that two of the protein concentrations approach the same steady state value. It is difficult to see a justification for that in the equations. Is it possible your spreadsheet does not correctly implement the differential equations shown?

Beyond that, I suggest that if Mathematica results (using adaptive time steps) differ from those in Excel using a fixed time step, then Mathematica is more trust worthy.

Attachments:

677785769-DK.nb

You are missing a semicolon on your first line between B2=1 and B3=1

If I fix that and then

sol=NDSolve[...yoursystem...];
Plot[Evaluate[{x[t],y[t],z[t]}/.sol],{t,0,100}]

Then I see your three plots. But the plots do not look very close to the plots you show from Excel. I am guessing that this may be because of the values for some of your coefficients when you tried to translate from discrete steps in Excel to derivatives in Mathematica. Perhaps you could estimate the derivatives by looking at the change in values from 0.0 to 0.000025 and the slope of those changes might be close enough to the derivatives that you need to give NDSolve

Bill,

I have attached a truncated version of my Excel file (the complete original file is 129 MB with 1e6 rows in it). You can see the formulas in the sheet that I used to simulate the concentrations of proteins x, y, and z from the differential equations for x, y, and z formation with respect to time. x, y, and z are the proteins made in the negative feedback system from the three genes, x, y, and z. All the parameters are defined in the Excel sheet too that are used in the three differential equations. Thank you for any help you can render to me.

Scott

Attachments:

Oscillating equa...xlsx

x,y, and z are the three different proteins that are made by the three-gene negative feedback loop and described by the three differential equations. Thanks for any help you can give me.