Line through the origin mod $1$ visits every sub-cube in $\mathbb{R}^n$.

So my understanding is that we are considering the line (segments) given by $$ \left\{ \begin{gathered} x_{\,1} = a_{\,1} \;t\;\bmod 1 = \left\{ {a_{\,1} \;t} \right\} \hfill \\ x_{\,2} = a_{\,2} \;t\;\bmod 1 = \left\{ {a_{\,2} \;t} \right\} \hfill \\ \quad \vdots \hfill \\ x_{\,n} = a_{\,n} \;t\;\bmod 1 = \left\{ {a_{\,n} \;t} \right\} \hfill \\ \end{gathered} \right.\quad \Rightarrow \quad \mathbf{x} = \left\{ {\mathbf{a}\;t} \right\}\quad \left| \begin{gathered} \;a_{\,k} \in \;\;\mathbb{Z}\,_ + \; \hfill \\ \;\gcd (a_{\,1} ,\, \cdots ,a_{\,n} ) = 1 \hfill \\ \end{gathered} \right. $$ where $\left\{x \right\}$ is the fractional part of $x$, i.e. $$ \left\{ x \right\} = x - \left\lfloor x \right\rfloor $$ and either it and the floor function are assumed to apply to vectors component-wise.

That premised, let us invert the scheme:
consider the space partitioned into a lattice of hyper-cubes with unitary sides and let the line "free" to run across the space, i.e. across the cubes of the n-D lattice.
We will individuate each cube according to the floor of the coordinates of the points it contains, that is according to $\left\lfloor \mathbf{x} \right\rfloor $, which are the coordinates of its lower vertex. That also means that each cube so individuated will include the faces with lower coordinates and exclude the others.

Starting from the origin, the line will traverse the cube $(0,\, \cdots ,\, 0)$ and continue across the others.
If in traversing a cube, the line segment within it has the same relative coordinates as the segment from the origin, then it will individuate a cycle, same as in the original scheme proposed.
Relative coordinates same as the segment from the origin, means a relative $(0,\, \cdots ,\, 0)$ point of origin of the segment, i.e. a point with integral coordinates.

Therefore we are looking for the values of the parameter $t$ such that $$ \left\{ {\mathbf{a}\;t} \right\} = \mathbf{0}\quad \Rightarrow \quad \left\{ {t\;\gcd (a_{\,1} ,\, \cdots ,a_{\,n} )} \right\} = 0\quad \Rightarrow \quad \left\{ t \right\} = 0\quad \Rightarrow \quad t \in \;\;\mathbb{Z}\; $$ We will have a repetition at each block $t \, (a_{\,1} ,\, \cdots ,a_{\,n} )$.
The number of blocks that the line traverses from the origin till the point at $t=1$ (excluded) will correspond, in the original scheme, to the number of distinct segments inside the cube.
We have a boundary cross whenever one of the components of $\left\{ {\mathbf{a}\;t} \right\}$ gets null, and $$ \begin{gathered} \left\{ {a_{\,k} \;t} \right\} = 0\quad \left| {\;0 \leqslant t < 1} \right.\quad \Rightarrow \quad \hfill \\ \Rightarrow \quad t \in \left\{ {0,\;\;\frac{1} {{a_{\,1} }},\; \cdots ,\;\frac{{a_{\,1} - 1}} {{a_{\,1} }},\;\;\frac{1} {{a_{\,2} }},\; \cdots ,\;\frac{{a_{\,2} - 1}} {{a_{\,2} }},\; \cdots \; \cdots \;,\frac{1} {{a_{\,n} }},\; \cdots ,\;\frac{{a_{\,n} - 1}} {{a_{\,n} }}} \right\} = \hfill \\ = \left\{ {t_{\,0} = 0,\;\;t_{\,1} = \frac{1} {{a_{\,n} }},\;\;t_{\,2} ,\; \cdots \;,\;\;t_{\,q} } \right\}\quad \left| {\,a_{\,n} = \max \left( {a_{\,1} ,\; \cdots ,\;a_{\,n} } \right)} \right. \hfill \\ \end{gathered} $$ We may assume, WLOG, that the components of $\mathbf a$ are in non-decreasing order. Since it might be that $\gcd (a_{\,k} ,a_{\,j} \ne 1$, for $k, j$ ranging from $1$ to $n$, some of the values for $t$ might be coincident, and the above expression shall be understood in fact as a set.
If $\gcd (a_{\,k} ,a_{\,j} = m$ then the corresponding sub-sets will have $m-1$ elements in common, so the cardinality of the set depends from the mutual $\gcd$s between the components of $\mathbf a$.

However we have that
number of crossings = number of different segments in the original scheme (apart the $0$) $q \leqslant -n+ a_{\,1} + a_{\,2} + \, \cdots + a_{\,n} $

If the components of $\mathbf a$ were in the arithmetic sequence $(b+1,\; \cdots ,\;b+n)$, then we could apply to the Farey sequence and tell that the the No. of values of $t$ be $$ q = \sum\limits_{b + 1\, \leqslant \,k\, \leqslant \,b + n} {\varphi (k)} $$ with $\varphi$ being the Euler Totient function.
In any case, the Farey's sequence tells us that the $t_{k}$ values , apart $0$, are symmetric with respect to $1/2$, and that $q$ is odd or even depending on whether $1/2$ is included or not in the set.

The example sketched above ($\mathbf a =(2,3,5)$) helps to visualize the process involved.
The line enters a cube of the lattice at $t_{k} \mathbf a$ and leaves it at $t_{k+1} \mathbf a$.
The relevant segment is translated back by $\left\lfloor {t_{\,k} \,\mathbf{a}} \right\rfloor $ to the reference cube at the origin. During the translation, the segment remains parallel to $\mathbf a$ and thus orthogonal to any of the vectors normal to $\mathbf a$.
The configuration inside the refence cube of the $q+1$ segments (including the last $t_{q},1$) will be symmetric with respect to the $(0,0, \cdots ,0)$ and $(1,1, \cdots, 1)$ vertices. They will lay on planes parallel to $\mathbf a$.
Let call $\mathbf b$ one of the vectors normal to $\mathbf a$, then

$$ \begin{gathered} \mathbf{b} \cdot \mathbf{a} = 0\quad \Rightarrow \quad \mathbf{b} \cdot \left( {t\,\mathbf{a}} \right) = 0\quad \Rightarrow \quad \mathbf{b} \cdot \left( {t_{\,k} \,\mathbf{a}} \right) = 0\quad \Rightarrow \hfill \\ \Rightarrow \quad \mathbf{b} \cdot \left\{ {t_{\,k} \,\mathbf{a}} \right\} = - \,\mathbf{b} \cdot \left\lfloor {t_{\,k} \,\mathbf{a}} \right\rfloor \hfill \\ \end{gathered} $$

and all a cascade of relations.
For our scope we shall seek for the vector $\mathbf b$ that minimizes the set of values corresponding to $\mathbf{b} \cdot \left\lfloor {t_{\,k} \,\mathbf{a}} \right\rfloor$ at the varing of $k$. This will minimize the number of planes containing the points, and thus provide the largest separation between them, which in turn will determine the maximum $M$ as requested.
Since the $t_k$ are rational, we may restrict the search to the vectors with integral components.
By conducting this exercise on a few examples, it comes evident and it is intuitive ( without pretending to constitute a rigorous proof ) that the vector in question is the one that

gives $\mathbf{b} \cdot \mathbf{a} = 0$ and has minimum modulus

and this is the vector provided by the Multidimensional Extended Euclidean Algorithm, which is connected to the Bezout Identity.
It is to remember in fact that in 2D the algorithm provides $\gcd (a,b) = n_{\,k} \,a + m_{\,k} \,b\quad 0 = n_{\,k + 1} \,a + m_{\,k + 1} \,b$ and either $\left( {n_{\,k} ,\;\,m_{\,k} } \right)$ and $\left( {n_{\,k+1} ,\;\,m_{\,k+1} } \right)$ have minimum modulus.
Finally $c= \mathbf{b} \times \mathbf{a}$ will give the separation among the segments in the same plane.