How to write this Matlab code in Mathematica

Question

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)));
end

Can anyone tell me how I can code this part?

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:

countOpenGates[randmat]

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], 
             Dimensions[mat]])]

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

 ListLinePlot[opengatepath]

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]