I have a small hiccup here. I am trying to write my Matlab code for a dynamical system in Mathematica. I have managed to get most of it except one part.

The Mathematica code for the system is

ctot = 2;
c1 = 0.185;
v1 = 80;
v2 = 0.11 ;
v3 = 27 ;
k3 = 0.146 ;
a2 = 0.2 ;
Gates = 500000;
OpenG = 50000;
alpha = 0.27;
k4 = 1.1 ;
v4 = 6 ;
Ir = 1 ;
d1 = 0.130 ;
d2 = 1.049 ;
d3 = 0.943 ;
d5 = 0.082 ;

sol = NDSolve[{c'[
     t] == (v1 (pp[t]/(pp[t] + d1))^3 (c[t]/(c[t] + d5))^3 (OpenG/
          Gates)^3) (ctot - (c1 + 1)*c[t]) - 
     v3*c[t]^2/(k3^2 + c[t]^2), 
   pp'[t] == v4 (c[t] + (1 - alpha) k4/(c[t] + k4)) - Ir*pp[t], 
   c[0] == 0.1, pp[0] == 2}, {c, pp}, {t, 0, 100}]

Plot[c[t] /. sol, {t, 0, 100}]

It seems all okay except for the part where I need to calculate OpenG

For the purpose of simplicity and to verify my Mathematica code I have assigned OpenG 50000. But this should be a varying parameter. In Matlab, this part of the code is to obtain the number of open gates(for example) driven by random conditions.

% randomly generate number of open gates
ini_states = unidrnd(2,100,3) - 1;
% a matrix of randomly generated zeros and ones. row corresponds
% to each gate, while column corresponds to the states in each gate's subunits.

alpha_h = 0.2*1.049*(2 + 0.13)/(2 + 0.943)*0.01; % a constant

openG = [length(nonzeros(all(ini_states,2))); zeros(100,1)];
% initial number of open gates. It checks for the rows that have all ones.
% if a row has colums with all ones then the gate is in the open state,

for i=1:100
    y = rand(100,3);  %probability that each gate's subunit will change its state

    beta_h = 0.2*0.1*c[t]; %replace c[t] with rand just for the purpose to checking the code in MatLab

    indces = (ini_states & y<beta_h) | (~ini_states & y<alpha_h);
    % 'indces' matrix  holds the position of the states that have been changed,
    %indicated by a logical 1. indexing the values of 'indces' to ini_states
    %and then inverting the values will give us the new state matrix.

    ini_states(indces) = ~ini_states(indces);

    %obtaining the new count on the number open gates
    openG(i+1) = length(nonzeros(all(ini_states,2)));

Can anyone tell me how I can code this part?

  • 3
    $\begingroup$ RandomInteger gets you a random integer. Select[#!=0&] selects nonzero elements of a list. ConstantArray[0, {dims}] creates a list whose elements are zero. $\endgroup$ Commented Jan 13, 2016 at 9:21
  • $\begingroup$ @PatrickStevens Thank you very much. I will look in those keywords. $\endgroup$
    – nashynash
    Commented Jan 14, 2016 at 0:22

1 Answer 1


This isn't a complete answer, just something to get you started.

I think it would useful to think of a function, say countOpenGates, which would look something like this:

countOpenGates[mat : {{(0 | 1) ..} ..}] /; MatrixQ[mat] := 
  Total[Times @@@ mat]

If "open" means "all subunits are open", i.e. have value 1, then their product is 1. And if any subunit is closed (value 0), then the product of the row is 0. Taking the product (Times) of each row by Applying Times to each row, that is, at level 1, can be done using the @@@ shorthand notation. Then to know how many gates are open, just Total the resulting list of ones and zeros. If you want to know the fraction of gates that are open, then you would instead write Total[Times @@@ mat]/Length[mat].

I should explain the syntax of the function definition, that is the bit on the left hand side of the := (SetDelayed).

The pattern (:) required is a list of lists (note the use of the .. for a repeated pattern) where each element is either a one or a zero. A single | means Alternatives -- it's different from Or (||) and is only relevant for pattern-matching. Since pattern-matching is one of the things that distinguishes Mathematica from Matlab, I recommend you have a look at the guide on patterns and the fast introduction for programmers on the topic.

You can set up your random gates like this:

 randmat = RandomInteger[1, {10000, 5}];

And count the number of open gates like this:


To flip the state of some subunits, you could do something like this:

flipstate[mat : {{(0 | 1) ..} ..}, prob_?NumericQ] /; 
  MatrixQ[mat] && 0 <= prob <= 1 := 
  Abs[mat - (RandomVariate[BernoulliDistribution[prob], 

The patterns to match mat are the same as before. Here, I've specified the probability of a subunit flipping as a parameter prob. Since that's just a BernoulliDistribution, you just need a random matrix with each element having the appropriate distribution. I then subtract that matrix from the original matrix, so that a flipped subunit is -1 if it was zero and 0 if it was one. Taking the absolute value Abs converts those -1s to 1s without messing anything else. (Note that Abs is Listable so you can get element-by-element absolute values by applying the Abs function to the whole list. No loops or special syntax required!)

The last line suggests that you want to flip the gates repeatedly and count how many are open at each iteration. This is an obvious job for NestList.

opengatepath = Last /@ NestList[
With[{flipped = flipstate[First@#, 0.3]}, 
  {flipped, countOpenGates[flipped]}] &, 
   {randmat, countOpenGates[randmat]},  200]

This gives just the countOpenGates result for each iteration, because I Map the function Last onto a result that contains both the gates and the count. You could of course have some little function that makes the number of iterations and the flip probability into parameters.

You can then do things like


Edit in response to OP comment

To solve the ODE for each value of opengatepath as calculated above, just do this:

Table[sol2[i] = 
 NDSolve[{c'[t] == (v1 (pp[t]/(pp[t] + d1))^3 (c[t]/(c[t] + d5))^3 (opengatepath[[i]]/
      Gates)^3) (ctot - (c1 + 1)*c[t]) -  v3*c[t]^2/(k3^2 + c[t]^2), 
 pp'[t] == v4 (c[t] + (1 - alpha) k4/(c[t] + k4)) - Ir*pp[t], 
 c[0] == 0.1, pp[0] == 2}, {c, pp}, {t, 0, 100}] , {i, Length[opengatepath]} ]

You can then do things like

Grid@Partition[Table[Plot[c[t] /. sol2[i], {t, 0, 100}], {i, Length[sol2]}], 5]
  • $\begingroup$ Thank you very much. This is helping me. But I wonder how can I pass the obtained OpenG value to the ODE in NDSolve. Also note that OpenG has to be calculated at every iteration of NDSolve. For example NDSolve is set from 0-100, so for the ith iteration we have to calculated the ith OpenG value and accordingly use it in the first ODE. $\endgroup$
    – nashynash
    Commented Jan 14, 2016 at 2:08
  • $\begingroup$ If the result of the ODE does not affect the next iteration of OpenG, then you can just use a Table that takes the i-th part of opengatepath. I don't actually understand what you are trying to do. t is continuous time in an ODE; it isn't discrete iterations. $\endgroup$
    – Verbeia
    Commented Jan 14, 2016 at 2:20
  • $\begingroup$ I've added something that I think is what you need, but I might have misunderstood you. $\endgroup$
    – Verbeia
    Commented Jan 14, 2016 at 2:27
  • $\begingroup$ The result of the ODE does affect the next iteration of OpenG. beta_h = 0.2*0.1*rand. I used rand here to make it easier to understand but this is suppose to be c[t]. So beta_h depends on c[t]. Sorry I should have included that as well in the first place. $\endgroup$
    – nashynash
    Commented Jan 14, 2016 at 2:32
  • $\begingroup$ Ok, this makes it more complicated, but essentially what you need to do is use NestList with a more complex body, including the ODE. I don't have time to do this for you today, but what I have above should get you started. $\endgroup$
    – Verbeia
    Commented Jan 14, 2016 at 22:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.