**2. Model and Simulation Details**

The CNTs are modelled as the ellipsoids of revolution with high aspect ratio, their semi axes are *b*<sup>1</sup> = *b*<sup>2</sup> *b*3. For the calculation the similar ellipsoids with semi axes *b*<sup>1</sup> = 2.5 nm and *b*<sup>3</sup> = 37.5 nm were taken. The representative unit cell with the volume of (*nb*3)3, where *n* > 1, is used.

The particles are introduced one by one into the unit cell by randomly generating Cartesian coordinates *X* for the centre of the ellipsoid and two spherical angular coordinates *θ* and *φ* for the orientation of the longest axis of the ellipsoid. Each time the non-overlap condition is checked, and the new ellipsoid is stored only if it does not intersect the walls of the unit cell and does not penetrate into any of the already existing particles. After the predefined number of particle is reached, the connectivity is checked along the selected direction using the Dijkstra's protocol [26].

The algorithm for the percolation computation is organised as a Tabu search method and it stops when the percolation is achieved [27–29]. Until the percolation is reached, the number of particles increases in each cycle. The detailed modelling procedure is described in our previous paper [30]. The periodic boundary conditions are applied for the calculation of the percolation concentration and conductivity.

The connection criteria for the nanotubes is obtained as follows. The dependence of the tunnelling resistance on the distance between the nanotubes given as [8,31]:

$$\rho = \frac{h^2}{\epsilon^2 \sqrt{2mb}} \exp\left(\frac{4\pi d}{h} \sqrt{2mb}\right) \tag{1}$$

where *b* is the potential barrier, and *d* is the distance, *e* and *m* denote the elementary change and electron mass, respectively. The reciprocal value 1/*ρ* is the tunnelling conductivity. The conductivity drastically decreases with the distance increment and reaches already small values of 10−<sup>4</sup> S/m for 2 nm separation. Thus, for the percolation computations the separation of 2 nm is considered as the connection criteria, while for the piezoresistivity computation direct values obtained from Equation (1) were used. In the last case, the tunnelling barrier of *b* = 0.75 eV.

To introduce anisotropy, the composite was considered as mechanically stretched. The probability density function (PDF) for the CNT angular distribution was derived using the following assumptions: (i) the mechanical deformation oriented along the *z*-axis, (ii) the centre of the ellipsoid keeps its position after the deformation, and (iii) the Poison's shrinkage in the perpendicular directions was neglected. The last assumption is justified, since the independence of the percolation value on the unit cell volume was demonstrated previously [30].

Thus, after the deformation, *X* and *ϕ* remain unchanged. The PDF for the *θ* angle was obtained as the modified function for the mechanically deformed composite filled with cylinders [32]:

$$\psi(\theta) = \frac{k^2}{[k^2 \cos^2 \theta + \sin^2 \theta]^{3/2}} \tag{2}$$

where *k* is the deformation coefficient introduced as the ratio of the final and initial lengths of the unit cell, *k* = *<sup>l</sup> l*0 . For *k* = 1, the PDF function (2) will return uniform distribution, while for *k* > 1 some non uniformity will appear (see Figure 1). In order to apply the function to spherical coordinates, the final distribution will be given taking into account the Jacobian as *ψ*(*θ*)*sin*(*θ*).

**Figure 1.** Probability density function (PDF) for the zenith angle for different deformations (*k*) for the carbon nanotubes (CNTs).

The conductivities are computed according to the following protocol. Firstly, the resistance of the nanotubes is taken as infinite, so the total conductivity of the composite is governed by the inter-tube tunnelling (Equation (1)). Next, the direction for the conductivity computation is selected and the nanotubes near the initial and final boarders are collected. The Dijkstra algorithm is used to trace the paths of the minimal resistance between the initial and final tubes. The array *R*(*t*, *l*) (where *t* and *l* stand for the initial and final indexes of the tube, respectively) is computed. To implement the periodic boundary conditions the array *B*(*l*, *t*) of the boundary resistances was introduced. The total resistance for the selected fixed *l* is follows:

$$Tot(l) = \left(\sum\_{t} \frac{1}{R(t, l)}\right)^{-1} + \left(\sum\_{t'} \frac{1}{B(l, t')}\right)^{-1} \tag{3}$$

And finally the total conductivity of the composite is computed as

$$\sigma = \left(\sum\_{I} \frac{1}{Tot(I)}\right) \tag{4}$$

The algorithm is implemented using PGI CUDA FORTRAN standards [33]. The big enough number of the realizations (500–600) was collected for each particular case.

## **3. Results and Discussion**
