Catalyst.jl - @observables usage

Specific Domains Modelling & Simulations
1

Hi, I have a question about using @observables inside @reaction_network -
How to write the observable contains both parameters and species? For example, what’s the correct syntax for writing Tree?

 function model()
   ....
   @species begin
         species_A1(t)
         species_A2(t)
   end

   @parameters begin
         p1
         p2
   end
   
   rn = @reaction_network begin

        @observables begin 
            total_A ~ species_A1 + species_A2 # No issue
            Tree ~ species_A1 * p1 - species_A2 * (p2-p1) # throw error
        end
         .....
   end

Thanks so much!

2

The @reaction_network macro is intended to have the complete model defined within it, so you need to put the @species and @parameters declarations within it. Otherwise they are completely decoupled from it (i.e. those species and parameters will not actually be seen in the macro and it will instead infer them if possible – however, it can’t infer unknown symbols within observables and that is why you get an error). The standalone @species and @parameters are intended to be used via the symbolic interface, i.e. where you separately create species, parameters, equations, and reactions one by one in code, and then manually pass them into a ReactionSystem.

So your example using @reaction_network should be:

 function model()
   rn = @reaction_network begin
       @species begin
             species_A1(t)
             species_A2(t)
       end

       @parameters begin
             p1
             p2
       end

        @observables begin 
            total_A ~ species_A1 + species_A2 # No issue
            Tree ~ species_A1 * p1 - species_A2 * (p2-p1) # throw error
        end
   end
   return rn
end
3

Hi Isaacsas,

Thank you for getting back to me. Please see the full function below. The reason I have them outside of @reaction_network is that if I move @species and @parameters inside @reaction_network, I need to move the two rates inside too. But then I got an error: @reaction_network input contain 2 malformed lines (the two line describing rate).

If this is the case, what’s the best practice please? Many thanks!

 function model()
   t = default_t() # i have to add this, otherwise it throws error
   
   @species begin
         species_A1(t)
         species_A2(t)
   end

   @parameters begin
         p1
         p2
   end

   A_rate = p1 * species_A1
   B_rate = p2 * p1 * species_A1
   
   rn = @reaction_network begin

        @combinatoric_ratelaws false

        @observables begin 
            total_A ~ species_A1 + species_A2 # No issue
            Tree ~ species_A1 * p1 - species_A2 * (p2-p1) # throw error
        end
        
       $B_rate, B => A
       $A_rate, 0 => A
       
   end
   return rn
end
4

You could make a function that calculates that rate and call it like:

rate(p,S) = p*S
rn = @reaction_network begin
            @species begin
                  species_A1(t)
                  species_A2(t)
                  A(t)
                  B(t)
            end

            @parameters begin
                  p1
                  p2
            end

             @observables begin 
                 total_A ~ species_A1 + species_A2 # No issue
                 Tree ~ species_A1 * p1 - species_A2 * (p2-p1) # throw error
             end
             
             rate(p1, species_A1), B => A
             rate(p1*p2, species_A1), 0 => A
        end

otherwise you have to use interpolation to make the macro aware of externally defined parameters/species like (note if a species only appears in an observable it won’t be added as a species – I assume in your real model you want to use species_A2 in some rate law or reaction so I’ve modified your code slightly)

function model()
          t = default_t() # i have to add this, otherwise it throws error
          
          @species begin
                species_A1(t)
                species_A2(t)
          end

          @parameters begin
                p1
                p2
          end

          A_rate = p1 * species_A1
          B_rate = p2 * p1 * species_A2
          
          rn = @reaction_network begin

               @combinatoric_ratelaws false

               @observables begin 
                   total_A ~ $species_A1 + $species_A2 # No issue
                   Tree ~ $species_A1 * $p1 - $species_A2 * ($p2-$p1) # throw error
               end
               
              $B_rate, B => A
              $A_rate, 0 => A
              
          end
          return rn
       end
1 Like