MVAR Simulation for Microbiota and Metabolites

MVAR Process

The MVAR model extends the traditional VAR model used to fit abundances over time, by incorporating metabolites. The model can be expressed as follows:

\[X_t = A \cdot X_{t-1} + \epsilon_{X,t}\]

\[S_t = B \cdot X_{t-1} + \epsilon_{S,t}\]

where:

• \(X_t\) is the vector of microbiota abundances at time \(t\),

• \(S_t\) is the vector of metabolite abundances at time \(t\),

• \(A\) and \(B\) are matrices representing the interactions within microbiota and between microbiota and metabolites, respectively,

• \(\epsilon_{X,t}\) and \(\epsilon_{S,t}\) are vectors of error terms, assumed to be normally distributed with mean 0 and a specific standard deviation.

Objective

The objective of this notebook is to simulate the dynamics between microbiota and metabolites using the MVAR model. Through simulation, we aim to demonstrate how specific interactions can be identified and quantified, providing insights into the causal relationships within these systems.

Simulation with VARsim.py for sVAR Models

The VARsim.py script has been adapted to support MVAR simulation, allowing us to specify metabolite interaction matrices and simulate the dynamics of microbiota and metabolites over time. This section demonstrates how to use the script for simulating and visualizing the interactions captured by our MVAR model.

Example Usage of VARsim.py for MVAR Simulation

The implementation of the VAR model in VARsim.py allows for detailed simulation of the dynamic interactions between microbiota and their produced metabolites. The script is designed to facilitate straightforward simulation of these complex biological systems. The following code examples provide a guide on how to leverage the VARsim.py script to simulate these interactions effectively:

# Importing the libraries
from mimic.model_simulate.sim_VAR import *
from mimic.model_infer.infer_VAR_bayes import *

WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.

First, let’s first generate synthetic data for microbiota using the sim_VAR class. We will then simulate the sparse interaction with a given matrix

# Simulate a VAR model

simulator = sim_VAR()

We can now use the following parameters to generate the metabolic data:

```json { “coefficientsM”: [ [ 0.0, -0.5, 0.0 ], [ 0.1, 0.1, -0.1 ], [ -0.2, 0.1, 0.3 ] ], “initial_valuesM”: [ [ 2 ], [ 0 ], [ 0 ] ] }

simulator.read_parameters(r'parametersS.json')
simulator.print_parameters()

simulator.simulate('MVARsim')

data = simulator.data
print(data)

[[ 1.00000000e+00  2.00000000e+00  0.00000000e+00]
 [-2.74161993e+00  1.86811048e+00  9.39518955e-01]
 [-3.87544515e+00 -6.57954979e-01 -2.16334993e-01]
 [-4.57029479e+00 -7.24143855e-01 -2.22807962e+00]
 [-4.88519131e+00  5.61638642e-01 -2.00909298e+00]
 [-5.47635413e+00  2.27072968e+00 -3.38864564e+00]
 [-6.74251671e+00  7.44873096e-01 -3.43832484e+00]
 [-8.60768672e+00  3.82065959e+00 -2.80719684e+00]
 [-9.36604578e+00  3.47894584e+00 -5.06874169e+00]
 [-1.00013258e+01  3.96982781e+00 -7.21746706e+00]
 [-1.06438315e+01  7.05567873e+00 -4.63275738e+00]
 [-1.23546385e+01  5.59630984e+00 -4.48170642e+00]
 [-1.09923513e+01  3.89930726e+00 -4.24766812e+00]
 [-1.24157574e+01  2.83688134e+00 -3.60737571e+00]
 [-1.27878037e+01  8.40982423e-01 -2.33450717e+00]
 [-1.00426173e+01 -5.44655146e-01 -3.16483959e+00]
 [-8.01431316e+00  2.36703403e+00 -3.28244696e+00]
 [-6.63400199e+00  4.44366070e+00 -2.15491020e+00]
 [-5.77061264e+00  1.64892443e+00 -1.64449167e+00]
 [-8.45696391e+00  1.82216559e+00 -3.51180747e+00]
 [-8.01972123e+00  3.23649646e+00 -2.80390688e+00]
 [-1.02138737e+01  6.58194510e-01 -2.70335987e+00]
 [-9.08694475e+00  6.02318281e-01 -2.73878937e+00]
 [-7.34899973e+00 -6.91122789e-01 -2.46676053e+00]
 [-5.09694871e+00 -1.79029896e+00 -2.91559430e+00]
 [-3.38511951e+00 -6.24026896e-01 -2.64341719e+00]
 [-2.31563951e+00  1.63621392e+00 -6.40097031e-01]
 [-2.30380257e+00  6.19049430e-01 -1.11274792e+00]
 [-3.38343865e+00  4.86553494e-02 -1.24290930e+00]
 [-3.18075309e+00 -1.19071554e+00 -2.03524968e+00]
 [-4.14619831e+00  1.08682219e+00 -1.50655710e+00]
 [-1.49220530e+00  4.18451374e-01 -1.87388195e+00]
 [-2.93283276e+00  4.10210417e-02  1.23018332e+00]
 [-1.61872411e-01 -1.70699796e+00 -2.41753619e+00]
 [-9.63421966e-01  1.77772933e+00 -3.63334807e-01]
 [-3.96364593e-01 -8.13755989e-01 -2.98042395e+00]
 [-2.44820198e+00  3.57536668e+00 -2.26514196e-01]
 [-4.06881334e+00  2.03228413e+00  4.18171616e-01]
 [-3.62360701e+00  5.67372520e-03 -2.15716865e+00]
 [-1.38169514e+00  8.46321789e-01 -2.74370175e-01]
 [ 3.45235322e-01  1.78630003e+00 -1.87148199e+00]
 [-7.17412231e-01  3.04534227e+00 -4.03323366e-01]
 [-1.12779346e+00  1.53832507e+00 -2.24561622e+00]
 [-2.75580218e-01  2.16015312e+00  8.63285740e-01]
 [-9.75696318e-01 -6.52515313e-01  7.98145013e-01]
 [-2.73434230e-02 -2.58510590e+00 -3.56623310e-01]
 [-6.93359529e-01 -7.99268697e-01  1.07855179e+00]
 [-1.51635301e+00 -2.34595637e+00 -9.33347578e-01]
 [-2.00114545e+00  8.76418566e-01 -8.07926292e-01]
 [-1.79893927e+00  1.70200852e+00  1.49594588e-01]
 [-9.44198194e-01 -2.22879736e-01 -1.10475418e+00]
 [ 9.12688273e-01  1.37597686e+00 -4.47951529e-01]
 [-4.34241000e-01  1.83759585e+00  1.20051236e+00]
 [ 8.98155844e-01 -5.78053871e-01 -1.20129535e+00]
 [-1.21416621e+00  1.90883823e+00  1.80390600e+00]
 [-6.65911892e-01 -1.12224694e+00  9.55910467e-01]
 [ 1.31355093e+00 -5.70717860e-01 -2.80991597e-01]
 [ 8.74414389e-01 -2.00004640e-01 -6.62668874e-02]
 [-1.93647083e-01  2.78344309e-01  9.71332354e-01]
 [-8.75243793e-01 -1.89406617e+00 -1.27982163e+00]
 [ 1.94422997e+00  8.22097533e-01  2.90928187e-01]
 [ 2.33918527e+00  2.86725675e+00 -4.83894731e-01]
 [-9.23861764e-01  3.18033474e+00 -8.40267412e-01]
 [-1.20445667e+00  2.28850824e+00  5.68445435e-03]
 [-1.92849404e+00  2.34547777e-01 -2.33123163e+00]
 [-3.34586841e+00  2.34672747e+00 -2.01624475e+00]
 [-3.34788898e+00  2.52728347e+00 -1.51991069e+00]
 [-5.88621155e+00  2.11865751e+00 -1.94292010e+00]
 [-4.42259210e+00  1.05958226e+00 -1.27045204e+00]
 [-3.29093169e+00  3.84238060e-01 -2.28538353e-01]
 [-2.45133035e+00  1.15183165e+00  8.83362910e-01]
 [-3.31899396e+00 -1.79750754e-01 -6.98593424e-01]
 [-3.15809829e+00  1.13371069e-01  3.13157538e-01]
 [-1.84507478e+00 -1.59492122e+00 -7.48393442e-02]
 [-1.11018373e+00 -1.12935477e+00  1.22305111e+00]
 [ 5.25118935e-01 -1.00252899e+00  2.72474536e+00]
 [ 8.91041452e-01 -3.27387727e+00  5.54369918e-01]
 [ 2.54016205e+00 -2.95045207e-01  8.79663907e-01]
 [-5.91362425e-03  3.13657194e-01 -2.71350319e-01]
 [ 1.58784451e+00 -7.41245056e-03 -1.26467183e+00]
 [-1.63228379e-01  1.05112773e+00 -6.86763147e-01]
 [-7.67824157e-02 -8.51985904e-01  2.16711169e-01]
 [ 1.24532251e+00 -1.38832751e-02 -6.01840099e-01]
 [ 1.33056733e+00  8.74031684e-01  1.70527321e-01]
 [ 2.90009069e-01 -1.06492738e+00  6.49401035e-01]
 [ 2.44697912e+00  9.70868702e-01  1.00064082e+00]
 [ 3.24444004e+00  2.57636384e-01  1.71452801e+00]
 [ 2.11206845e+00 -3.62243060e+00  1.80819859e-01]
 [ 3.32046359e+00 -2.23478665e+00 -6.32554111e-01]
 [ 1.91771241e+00  5.92653006e-01  1.09094443e+00]
 [ 1.24196350e+00 -3.51099096e-02  1.31372532e+00]
 [ 3.41808456e+00 -7.95776256e-01  1.20485825e+00]
 [ 3.92122967e+00  5.33326905e-01  2.88426353e+00]
 [ 5.70629672e+00 -3.02120823e+00  2.04426193e+00]
 [ 4.34835802e+00 -3.44844280e+00  2.51866752e+00]
 [ 4.88795221e+00 -2.56953466e+00  3.40203007e+00]
 [ 2.46271608e+00 -3.00455979e+00  2.06221286e+00]]

dataM = simulator.dataM
print(dataM)

[[ 2.          0.          0.        ]
 [-1.71547582 -0.65481951  2.59624547]
 [-1.09281769  0.08538762  0.77389912]
 [-0.13615964 -1.5354976   1.24711877]
 [-0.14448161  0.56095261 -0.2806991 ]
 [ 0.92287517  1.27321583  2.04177705]
 [-0.45398762  1.1010254  -1.30612694]
 [ 1.11843095 -0.09460795  1.47248469]
 [-2.60287057  0.20180641  2.50611568]
 [-1.13735474  0.30000163  1.79083028]
 [-0.89950133 -0.51800366  0.64938282]
 [-1.89834604  1.39916576  0.34556553]
 [-2.36499338 -1.31294636  1.0109132 ]
 [-1.22786393  1.54011497  0.08665203]
 [ 0.21287421 -1.02482448  2.33154482]
 [-0.13694572 -0.72161602  2.67291959]
 [ 0.551344   -1.3180056   0.5106264 ]
 [-2.53286945  1.71807724 -0.43080879]
 [-5.06772035 -0.41141136  2.86371985]
 [-1.60679158 -1.38431238  0.66713068]
 [-2.54694948 -0.3046455   3.19383167]
 [ 0.18517222  0.06649388  2.50467695]
 [ 2.45229018 -1.4116436   0.5935728 ]
 [-2.01465394  0.01196858 -0.88094006]
 [ 1.65873796  0.0124109   1.69210477]
 [ 1.57350599  0.24715309 -0.1629008 ]
 [ 2.40788051 -2.07699465 -1.69999898]
 [-2.70443573  0.33045716 -0.8862661 ]
 [-1.41919759 -0.97054     1.08847541]
 [-0.0523808   0.90202768 -0.92258785]
 [-0.83882107 -3.20295881  1.32960931]
 [ 1.3717529  -1.12749731 -0.93421857]
 [ 0.12626437  1.40222212  1.1185425 ]
 [ 0.7644137  -0.68819848  0.63370505]
 [ 0.7491966   0.18650051 -0.37206307]
 [-0.22196378  1.17266511  1.68208129]
 [-0.41105076 -2.29716532 -2.24989462]
 [-0.69919084 -1.72311334  1.89470891]
 [-1.37326902  0.79219683  3.0524694 ]
 [-0.2725896   0.57453144 -0.67314609]
 [-2.28621089 -2.69702486 -1.22034672]
 [-0.74155141 -0.06310745  2.30530771]
 [-2.68180988 -0.2245736   0.4272798 ]
 [-0.67418902  0.36945024  1.66993208]
 [-2.22088253  1.11116639  2.38867991]
 [-1.66362559  0.64886725  3.96577563]
 [-0.52614694 -0.01561298 -1.27477531]
 [-0.67813177 -1.12767118  0.52201962]
 [ 0.7868318   0.01386993 -0.59812055]
 [-1.35872278  2.5470283   0.47540161]
 [-1.64839699  1.96842778  0.25308793]
 [-1.69145186 -1.01958321 -0.07022641]
 [-0.5253923  -0.3107555  -0.6357894 ]
 [ 0.32707384  0.61549412  1.47981172]
 [ 1.73945456 -0.11107721  1.15113383]
 [-0.61043744  0.838662    0.43314448]
 [ 1.22234661 -1.72915457 -0.6666005 ]
 [-0.34502077  0.50049242  1.29582681]
 [ 2.62302985 -0.74075048  0.66458648]
 [ 0.8550264  -0.17359211 -1.23330509]
 [-0.92376831  1.19705193  1.39213048]
 [ 2.2871204   1.03039487  1.63918216]
 [-3.93569978 -1.98814866  0.18299848]
 [-1.89531664  2.03271694  3.44616531]
 [-2.03805636 -0.83313728  3.30945167]
 [ 1.13856898 -0.5335198   1.11028191]
 [-1.39162683 -1.0730105  -0.44182141]
 [-1.01703595 -0.9312178  -0.11209957]
 [-3.61304913  0.51480224  0.41952081]
 [ 1.38162618 -1.00200162  2.10563739]
 [-0.56251345 -0.41258432  1.10731338]
 [-0.12139003 -0.69779232  0.14475198]
 [ 0.74617752 -1.20058043 -0.54507379]
 [-0.4578581  -0.74702911  2.46831357]
 [ 1.57326133 -0.64609146 -0.97262917]
 [ 1.93186568 -0.301137    0.75833113]
 [ 1.26838163 -0.84396163  0.87977525]
 [ 0.73962047 -1.31285863 -2.20205594]
 [-1.16726445  1.05190592  0.26519363]
 [ 1.29718239 -0.98127527 -1.05932124]
 [ 2.23148411 -1.64703832  0.99898002]
 [ 1.1725964  -1.18378475  0.62861743]
 [-0.16199045 -0.85987166  0.43060607]
 [-1.8706776   0.10712541 -0.02292976]
 [-0.28308915 -1.52885834 -2.17706103]
 [ 3.14620807 -0.46565088  0.86492883]
 [ 0.54059912 -0.95402088  0.13549907]
 [-2.72129915  1.47102116 -0.65130233]
 [ 1.86874065  0.52150719 -2.23845514]
 [-0.56391762 -0.98493444 -3.1367987 ]
 [ 0.14579885 -0.00604977 -1.08298192]
 [ 2.84303666  1.74024124 -1.85194237]
 [ 1.51520858  3.00194406 -2.09104052]
 [-0.71779993  1.09694692  0.18716896]
 [-0.06464283  1.16371818  0.53827689]
 [ 0.03486117  0.44199979 -2.78265761]
 [ 1.22797718  0.02860057 -0.8845344 ]]

infer = infer_VAR(data=data, dataS=dataM)

infer.run_inference_xs()

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [X0h, S0h, Ah, Bh, sigma]

100.00% [6000/6000 00:35<00:00 Sampling 2 chains, 0 divergences]

Sampling 2 chains for 1_000 tune and 2_000 draw iterations (2_000 + 4_000 draws total) took 42 seconds.
We recommend running at least 4 chains for robust computation of convergence diagnostics

Results saved as:
NetCDF file: model_posterior_xs.nc
Data file: data_xs.npz

infer.posterior_analysis()

# we can also plot the results with the original paremetr values that we used to simulate the data

A = np.array([[0.8, -0.2,  0.3],
              [0.3,  0.5, -1.],
              [0.2, -0.1,  0.4]])

B = np.array([[0.,  -0.5,  0.],
              [0.1, 0.1, -0.1],
              [-0.2, 0.1, 0.3]])

infer.posterior_analysis(A=A, B=B)

infer2 = infer_VAR(data=data, dataS=dataM)

infer2.run_inference_large_xs()

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [c2_A, tau_A, lam_A, Ah, c2_B, tau_B, lam_B, Bh, sigma]

100.00% [24000/24000 02:26<00:00 Sampling 4 chains, 2,446 divergences]

Sampling 4 chains for 2_000 tune and 4_000 draw iterations (8_000 + 16_000 draws total) took 159 seconds.
There were 2446 divergences after tuning. Increase `target_accept` or reparameterize.

Results saved as:
NetCDF file: model_posterior_large_xs.nc
Data file: data_large_xs.npz

infer2.posterior_analysis(A=A, B=B)