I have done some more research on how to simulate with three coupled differential equations, but the code i have written in the attached Mathematica notebook is still not working. I am open to any advice to successfully simulate these three differential equations in Mathematica.

Attachments:

coupled three di...nb

Bill,

Thank you immensely for the help on this! Much appreciated. I did have a few of the parameters not set correctly like K2=5e7 not 1e7 and the starting condition for x[0] should have been 1e-6 not 0. However, the bad news is that even after setting those parameters and initial conditions correctly and setting the range from 0 to 0.000025 the curves do not match the Excel curves at all as you pointed out.

I think your suggestion is a good one if I understand it correctly. Do you mean NDSolve in Mathematica can be given a time increment like 0.000025 to do the numerical solution with or are you suggesting something else? At this point, I am not sure I am capable of doing either of your two potential solutions because I am not an expert in Mathematica.

Could you provide a little more detail on any ideas you have of overcoming the differences seen in the curves between Mathematica and Excel? It is hard for me to believe Excel would be better than Mathematica at this calculation, and it is more likely my incompetence in Mathematica may be causing the differences we are seeing in the curves. I am pretty confident that Excel (while taking forever with one million rows to do the calculation) is producing the correct curves.

Scott

David,

Thank you for all the help. I will take a look at the notebook you attached. Your analysis is much appreciated. To answer your question about two of the proteins reaching the same steady-state--that is correct because two of the proteins (x and z) have all the same parameters so upon reaching equilibrium (steady-state in essence) their protein concentrations should be the same. For some reason, protein x and z are not reaching the same steady-state level in Mathematica. In regards to your second question, I did find an error in my Excel sheet (I attached a truncated version yesterday, which has the mistake in it). However, the mistake did not change the oscillatory nature of the curves and the correct curves are almost identical those you saw in my post 2 days ago, so I am still perplexed why I am not seeing the same behavior in Mathematica. The curves should oscillate and two of the proteins should reach the same steady state. I have attached a corrected truncated version of my Excel sheet in case you care to look at it. I can find no mistake in the way I did the analysis in Excel other than it is not eloquent and extremely slow due to 1 million rolls in Excel. If I understand what you did in your notebook, I believe you just told Mathematica how large of a time increment to take when solving the coupled differential equations; is that correct? If that is correct, it is even more puzzling to me because the increments appear very small, which is similar to what I used for time (0.000025 time increments) in Excel.

Scott

Attachments:

Oscillating equa...xlsx

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

David,

yes, I screwed up the evaluation of that notebook--our messages crossed on the wire. It will not be the first and last time I do that I am sure! Mathematica is a real insider's software in my opinion.

I have attached the pages you requested, see Fig. 8-80. You can forget the hill coefficient exponents in the text. Hill exponents are for systems with positive or negative cooperativity (where binding of one ligand either increases the affinity or decreases the affinity of the next ligand binding event). This current system is complicated enough without throwing Hill exponents into the mix. That is easy to do later once the foundation is laid, which I think we are finally there now from the help of users on in this community like you. I am very grateful for the help on this.

Scott

Attachments:

Molecular Biolog...pdf

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

Hi Scott, I’ll have a look at this tomorrow.

David

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

Do you have the original differential equations and parameters that you used to develop your discrete model? That might make it easier to produce the differential equations that MMA can accept, solve and plot. Or did you develop the discrete model from scratch?

Both Excel and MMA can be used to approximate things like this. The work is trying to translate the information into a form that each of them can use.

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.

Bill,

Thank you for all the work! I chose Kx=1e8, Ky=5e7, and Kz=1e8, which could be realistic affinities (1/K) of 10 nanomolar and 20 nanomolar. Figure 8-83 is for the two-gene protein negative feedback loop. i chose the three-gene negative feedback system because the more genes the more the system oscillates. A five gene loop would have been even better, but I did not want to get carried away. You can leave the Hill coefficients (hx or y) out, if you would like, because a system with no cooperativity has hx or y=1. This system is complicated enough without more complexity. I choose all those 12 parameters arbitrarily though trial and error to create an oscillation the best I could in Excel. It will be very informative to see what happens in Mathematica when you enter your equations.

Scott

Hi Bill, I think we're both having a lot of fun with this!

In your code above I believe there is a factor missing in each equation. For example, in the first equation (1+Kz) should be (1+Kz z[t]). The [z] in the spreadsheet graphic is a bit confusing. I think it's how chemists write "concentration of z."

Best, David

Thank you, Closer. Still not the same, but this looks almost the same as what David posted 21 hours ago

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

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.

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,

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

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 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

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.