How to plot some but not all of the solutions generated by ParametricNDSolve

My question is about the syntax of plotting solutions of ParametricNDSolve but not all. Specifically how can I plot just $x_1,x_2,x_3$ not $x_1,x_2,x_3,i$. My attempt is below(no error messages)

A = {{2, 5, 0.5}, {0.5, 1, mu}, {1, 0.5, 1}};
tmax = 100;
sol = ParametricNDSolve[{x1'[
t] == -(1 + x1[t]) (A[[1, 1]] x1[t] + A[[1, 2]] x2[t] +
A[[1, 3]] x3[t] + e i[t]),
x2'[t] == -(1 + x2[t]) (A[[2, 1]] x1[t] + A[[2, 2]] x2[t] +
A[[2, 3]] x3[t] + e i[t]),
x3'[t] == -(1 + x3[t]) (A[[3, 1]] x1[t] + A[[3, 2]] x2[t] +
A[[3, 3]] x3[t] + e i[t]),
i'[t] == x1[t] + x2[t] + x3[t] - i[t],
x1[0] == x10, x2[0] == x20, x3[0] == x30, i[0] == 1},
{x1[t], x2[t], x3[t], i[t]}, {t, tmax}, {x10, x20, x30, mu, e}];
sol
Manipulate[ParametricPlot3D[
Evaluate[{x1[x10, x20, x30, mu, e][t], x2[x10, x20, x30, mu, e][t],
x3[x10, x20, x30, mu, e][t]} /. sol],
{t, 0, tmax},
Ticks -> None, PlotRange -> All, AxesOrigin -> Automatic,
AxesLabel -> {X1label, X2label, X3label},
BoxStyle -> Dashing[{0.02, 0.02}] , PlotRangePadding -> None,
AxesEdge -> {{1, 1, 1}, None, Automatic},
AxesStyle -> Thickness[0.005]],
{{x10, 1, "x1"}, 0.01, 10},
{{x20, 1, "x2"}, 0.01, 10},
{{x30, 1, "x3"}, 0.01, 10},
{{mu, 71/48, "mu"}, 0, 3},
{{e, 0.5, "e"}, -1, 1}]


Notes

• This is in part based on code from a previous post. My question was in that previous post was answered. See the post for a previously working version for 3 odes.
• If parametricNDSolveValue can be used instead and solve my problem, help with the proper syntax would be appreciated.
• What is A? It seems it is a matrix. ParametricNDSolveValue can be used. Commented Jun 14, 2018 at 0:48
• You may try /. sol[[1;; 3]]] Commented Jun 14, 2018 at 0:48
• $A$ is a matrix. I used ParametricNDSolveValue previously (see link in post) put I could not get it to work. Replaced /.sol with /. sol[[1;;3]]] it did not work, I got the same output.
– AzJ
Commented Jun 14, 2018 at 0:54

You're problem was that you solve for x1[t], x2[t] etc instead of x1, x2 etc. Fixing that yields:

sol = ParametricNDSolve[
{
x1'[t]==-(1+x1[t]) (A[[1,1]] x1[t]+A[[1,2]] x2[t]+A[[1,3]] x3[t]+e i[t]),
x2'[t]==-(1+x2[t]) (A[[2,1]] x1[t]+A[[2,2]] x2[t]+A[[2,3]] x3[t]+e i[t]),
x3'[t]==-(1+x3[t]) (A[[3,1]] x1[t]+A[[3,2]] x2[t]+A[[3,3]] x3[t]+e i[t]),
i'[t]==x1[t]+x2[t]+x3[t]-i[t],
x1[0]==x10,
x2[0]==x20,
x3[0]==x30,
i[0]==1
},
{x1,x2,x3,i},
{t,tmax},
{x10,x20,x30,mu,e}
];
Manipulate[
ParametricPlot3D[
Evaluate[Through @* (Through@{x1,x2,x3}[x10, x20, x30, mu, e]) @ t /. sol],
{t,0,tmax},
Ticks->None,
PlotRange->All,
AxesOrigin->Automatic,
AxesLabel->{X1label,X2label,X3label},
BoxStyle->Dashing[{0.02,0.02}],